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SUMMARY 


The contents of this report can be divided into five basic parts. The first 
part presents the state of the art in shell of revolution analysis, and briefly 
presents the theoretical advantages and disadvantages of the main numerical methods. 
The second part (Chapters 1-3) describes the static analysis pertinent to the 
STARS-2S program, and, in addition, presents the necessary equations for extension 
into nonlinear analysis of unsymmetric loading. The third part (Chapter k) deals 
with the analysis of the classical buckling loads of shells of revolution under 
axisymmetric loadings ( STARS -2B program), and unsymmetric loadings (not programmed). 
The fourth part (Chapter 5) deals with vibrations. The vibration and critical speed 
analyses involving axisymmetric prestress are programmed in the STARS -2V program 
while the analyses involving unsymmetric prestress remain unprogrammed. The final 
part (Chapter 6) presents several analytical studies performed with the STARS 
programs, where certain formulation advantages, or discrepancies with other analyses 
were -uncovered. 

It will be noted that the present eigenvalue formulation for the axisymmetric 
problem has several advantages over previous formulations. Possible similar 
advantages for the unsymmetric problem are discussed in Chapter 4, however there are 
many unanswered questions in this more complex area. Several comparative studies 
need to be carried out before a more accurate picture of even potential advantages 


can be drawn. 
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INTRODUCTION TO NUMERICAL ANALYSIS OF SHELLS 


Recent innovations in digital computer technology have enabled 
designers to analyze shell structures of complex configurations without 
unduly restrictive approximations. A number of versatile computer 
programs, based on various methods of analysis are presently available 
for the analysis of shells of revolution. With the exception of closed form 
solutions, the most commonly employed numerical methods of analysis 
will be briefly discussed and compared in the sequel. 

Finite Difference Method : This method was first employed by Radkowski, 
et al, [1, 2]* in analyzing layered shells of revolution subjected to 
axisymmetric loading. In these references, the two second-order 
differential equations of the theory presented by Reissner [ 3] were solved 
by using central differences with a constant mesh. The resulting 
simultaneous algebraic equations were solved by Potter's method £4] . 

Reference 1 is one of the few references wherein the finite difference • 

method is employed in analyzing branched shells. Only"Y branches" are 
considered. Sepetoski, et al. [5] employed a similar approach to study problems 
which do not involve shell branching. However, the simultaneous equations 
resulting from the application of finite differences were solved by a 
non-matrix Gaussian elimination technique. Moreover, in this reference, 
the effects of the grid sizes and of the error accumulation in the Gaussian 
elimination technique on the convergence of the solution were studied. 

Budiansky and Radkowski [6] extended the technique presented in Reference 1 j 
to the analysis of shells of revolution (including shells with " Y branches ") 
subjected to unsymmetric loading, using a shell theory presented by Sanders [7], 


* Numbers in brackets refer to the bibliography at the end of the report. 
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The load was expanded in Fourier series and a variable finite difference 
grid pattern was used. Greenbaum [8] presented a refinement for the analysis 
of shells closed at the apex. Capelli, et al. [9j included the effect of 
shear deformation and extended the technique of Reference 6 to shells 
of variable thickness in the circumferential direction. The variable 
thickness was expanded in a Fourier series, and the thickness harmonic 
amplitudes were coupled with the load harmonic amplitude s,r necessitating 
a modification of the procedure of Reference 6 to provide coupled harmonic solutions, 
Techniques similar to Reference 6 were applied in References 
10, 11, 12 to nonlinear problems. Hubka [10] presented an analysis of 
orthotropic, axisymmetrically loaded, shells of revolution, wherein non- 
linear terms were included in the equilibrium equations. Problems in- 
volving shells with two boundaries, subjected to distributed and concentrated line 
loads, were considered using Reissner's theory. The two second-order 
differential equations were solved using finite difference approximations 
with a variable grid. Options for forward, central, and backward differences 
were included. The resulting simultaneous equations were solved by a 
Gaussian elimination method specialized to banded matrices [ 13] . Nonlinear 
problems were solved by first solving the linear problem, and employing 
the computed in-plane stress resultants as input for the solution of the 
nonlinear equations. A similar method was utilized by Wilson and Spier 
[14, 15] to solve axisymmetric nonlinear shell problems using the equations 
of Reissner [ 16] . Moreover, they modified Potter's technique by adding 
an iteration procedure for solving the simultaneous nonlinear equations. 

In this iteration process the matrices generated by the solution of the linear 
equations were used as the first approximation in the solution of the non- 
linear problem. The procedure was repeated using the output of the previous 
iteration. This method converges very slowly for problems with significant nonline a ] 
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effects and it is necessary to use small load increments for higher values 

of the load. It was recognized by the authors that a Newtonian -type iteration 
technique [ 17, 18] would be more suitable inasmuch as its convergence 
could be proven. Busbnell [19] utilized the Newton -R aphson method to 
solve nonlinear problems of orthotropic, eccentrically stiffened shells 
of revolution. The effect of the stiffeners was treated by-" smearing"* the 
properties of the closely spaced rings and stringers over their spacing. 

Critical loads at buckling. for a shell subjected to axisymmetric loads 
(for cases in which the buckling mode is also axisymmetric) are established 
with this technique. The buckling load in established as the last load for 
which the Newton iteration method converges. 

The same technique, with a more approximate iteration procedure 
was first applied to shells of revolution subjected to unsymmetric loading 
by Ball [20, 21]. Inasmuch as the loading is expanded in a Fourier series, 
the nonlinear terms included in the Sanders shell theory [7], used in Reference 
20, couple the Fourier harmonics. The equations are uncoupled by first 
solving the linear problem, and using the results to obtain and introduce 
into the equations a numerical value for every nonlinear term. The resulting 
uncoupled equations are solved by the methods of Reference 6, each solution 
providing the numerical values for the subsequent iteration. The buckling 
load is defined as the last load for which the solution converges. Stability 
cases where the buckled shape is out-of-phase with the applied load cannot 
be considered. With this technique, it should be noted that the number of 


* See Appendix A for a further discussion of this technique. 
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sets of uncoupled equations which must be solved is not constant, but 
increases as convergence progresses, inasmuch as coupling terms 
involving more harmonics are evaluated numerically as the analysis 
progresses. Greenbaum and Conroy [22] found that for certain problems 
the technique of Reference 20 did not converge, and hence utilized a 
Newton iteration procedure. Moreover, it was determined that the 
formulation of the problem in terms of sets of four second-order 
differential equation results in numerical inaccuracies in the finite 
difference solution. For this reason the formulation of the problem in 
terms of sets of eight first-order equations was employed. Complexities 
of this nature in the use of finite differences are discussed in Reference 23. 

The technique presented in Reference. 22 has certain inconsistencies. In 
using the Newton iteration procedure, the cross -coupling terms of the 
various harmonic amplitudes are dropped in order to work with uncoupled 
sets of equations. In addition, some terms from the Sanders' shell theory 
are also dropped. Finally, the analysis of Reference 22 cannot be used 
for post -buckling analysis or for establishing the critical load at buckling 
where the buckling mode is out-of-phase with the load. 

Eigenvalue problems of vibrations and stability of shells have also been 
solved using the finite difference method. Cooper [24] has investigated the 
stability and natural vibrations of shells of revolution with two end 
boundaries, under axisymmetric prestress, using the methods of Reference 
6. The nonlinear theory presented by Sanders [ 7] is used to establish 
the equations of both the prestress state, and the incrementally disturbed 
state. In the vibration problem the perturbation displacements are assumed 
proportional to e lW *» and the frequency determinant is evaluated for a 
sequence of assumed values of the frequency, until the correct frequency 
is obtained as the one for which the determinant vanishes. Central differences 



with a constant mesh are employed to reduce the four second -order 
differential equations to algebraic equations. The latter are solved by 
Potters* method. This method is modified [ 25] in order to avoid spurious 
changes of sign of the determinant. The numerical procedure for establishing 
the buckling load is identical to that for establishing the frequency in the 
vibration problem. In the stability analysis, the prebuckling stress -resultants 
and deformations are considered. In the vibration analysis, the rotational 
inertia is omitted. Rossettos and Tene [26, 27] applied a very similar 
technique to the analysis of layered and orthotropic shells. The sole difference 
■was that they utilized second -order finite difference approximations in order 
to consider boundary conditions at the two ends with an accuracy consistent with 
that used for the analysis of the remainder of the shell. Heard and Fulton [ 28] 
also used an almost identical technique with a variable finite difference mesh. 

In References 29-37, a program for solving shell stability problems is presented. 

This program includes the effects of eccentric reinforcement on shells 

(smearing technique is used). Both the prebuckled and stability equations 

are treated by applying central differences to two fourth-order differential 

equations, and solving the resulting set of algebraic equations by use of 

the method presented in Reference 38. If the prebuckling state is established 

on the basis of a linear analysis, the stability equations may be conveniently 

approximated by the form ([A] + X [B])jx} = 0. Thus, the power method 

[39] may be employed to establish the critical load at buckling, avoiding 

the uncertainties inherent in the determinant evaluation method. In solving stability 

problems with the aforementioned program (BOSOR) the prebuckled state 

may be established by the nonlinear analysis of Reference 19. In such 

an eventuality, the determinant evaluation method is utilized. 
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Two-dimensional finite difference solution techniques have been used 
in solving vibration and stability problems in References 40 and 41. 

These two-dimensional techniques result in very large matrices, and 
therefore should be employed only in the absence of an alternative method 
(e.g. to investigate a shell cutout problem). 

I 

Finite Element Method : In applying the finite element method which 
is actually an application of the Rayleigh-Ritz numerical technique, to 
the analysis of shells of revolution, two distinct types of elements have 
been employed: the discrete triangle or quadrilateral; and the revolved 
conical or curved elements. Each of the foregoing elements have special 
advantages and specific areas of application- Generally, in the 
application of finite elements to shell analysis, the following two basic 
questions remain to be resolved: 

1. What is the effect of the geometric approximations between the 

. : , I ■; . ; - ■ - ; 

elements and the actual shell surface, on the solution of the problem. 

2. Is it necessary to explicitly include rigid-body motion terms in the dis- 
placement function employed in the solution. 

The first conical frustum element was introduced by Meyef and Harmon [ 42] 
The deformation of this element included membrane and bending components, 
and continuity of slopes and displacements was enforced on the inter -element 
nodal circles. The deformation of the element was established on the basis 
of an analytical solution for a cone loaded solely along the edge, and shell 
problems involving edge loading only were solved using the force method. 

Grafton and Strome [ 43] derived the matrices for a conical element using 
the displacement method. Simple polynomial forms were utilized to 
represent the deformation state. Friedrich [ 44] used a thick conical 
element whose deformation included shear deflections established by sample 
beam theory. Another conical element using an analytical edge-Ibaded 
cone solution for the deformation^ of the element was formulated by 

6 



L.U, et al. [45], using the displacement method, Percy, et al. [46], and 
Wilson [47], wererfche first to apply conical elements to solve problems 
of shells subjected to un symmetric loads. The non-symmetry was 
analyzed by use of Fourier series. In Reference 46 the effects of 
utilizing various order polynomials to represent the displacement functions, 
were investigated. 

A thorough study of the conical frustum element was performed by 
Percy, et al. , [46] and by Jones and Strome [ 48], They encountered 
considerable disadvantages in the use of this element for the analysis of 
doubly-curved shells. For predominantly membrane problems, due to the 
change of angle between adjacent elements, the kinks at the nodes produce 
calculated meridional bending moments where none exist. For problems 
wherein the rotations at the boundary are unconstrained, erroneous 
displacements of the boundary are calculated,: Moreover, it is essential 
to use extremely short cone elements a considerable distance in from the 
end boundaries. As indicated in Reference 48, optimization of shells by 
careful variation of the thickness or geometry is not possible when conical 
frustum elements are used. In static problems using conical elements , 
sofrie of the 1 errors cancel each other. However, this cannot be anticipated 
in the case of dynamic analyses where rapid meridional membrane stress 
variation may occur. To overcome the foregoing difficulties, Jones and 
Strome [49] introduced a doubly-curved revolved element where the coordinates, 
slope, and hoop principal radius of curvative, were continuous 
functions throughout, and, at the nodal circles, were identical to those of 
the actual shell. However, the meridional radius of curvature was a 
discontinuous function at the nodal circles, A similar element with special 
refinements for application at the apex of a shell was developed by Stricklin, 
et. al. [50], Khoja steh-Bakht [ 51] developed the matrices for a doubly-curved 


7 



revolved element matching the tangents and curvatures of adjacent elements 
at the nodal circles. A geometrically ’Still more refined revolved element was 
produced by Brombolich and Gould [ 52] . 

Stricklin, et.al., used the elements of Reference 50 to solve 
nonlinear problems of axi symmetrically loaded shells whose thickness 
and material properties were both axisymmetric [ 53] and non-axisymmetric } 
[54], The non-axisymmetric, nonlinear problems resulted in a set of • 

simultaneous equations involving coupled harmonic amplitudes similar to 
that solved in Reference 9 (on the basis finite difference techniques): 

It has been shown that it was not necessary to include explicit rigid body 
displacement terms [ 50, 54, 55, 56] in the deformation functions of these revolved 
curved elements. 

Dynamic and stability problems were also solved by the use of 
these revolved elements. Klein and Sylvester [ 57] and Bacon and 
Bert [ 58] solved shell vibration problems by using conical frustum elements. 

The latter authors included transverse shear deformation and rotatory 
inertia, and analyzed sandwich shells. Mode shapes and frequencies of 
orthotropic shells of revolution were found by Leimbach, et.al. [ 59], and 
subsequently by A delman, et.al. [60, 61] using elements better fitting the 
geometry of the shell. The mode shapes of sandwich shells were calculated 
by Abel and Popov [62], and nonlinear vibration problems were considered 
by Stricklin et. al. [63] . Navaratna, et. al. [64] applied both the cdnical 
frustum and the curved revolved elements to the solution of linear buckling 
problems, for shells of revolution subjected to axisymmetric loads. Each 
possible Fourier harmonic mode of buckling must be checked to find the 
lowest eigenvalue. 

It should be noted that the principal advantage of utilizing revolved 
finite elements instead of finite differences is that arbitrarily branched 
shells of revolution may be analyzed in a routine manner by using revolved 
finite elements. 
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Shells were also analyzed using discrete triangular or quadrilateral 
elements. Flat triangular elements were developed for this purpose by 
Melosh [65]. These elements accounted for membrane and bending 
flexibility, and may be used to analyze shells with material properties ranging 
from isotropic to aelotropic. Zienkiewicz, et. al. [66] applied triangular 
plate elements to the analysis of thin arch dams, notwithstanding the continuous 
interest in these applications [67], flat elements may be more inadequate 
than the conical frustra [48] in the representation of curved surface^ 
particularly in predominantly membrane stress areas. Solutions obtained 
using a number of flat elements were tested for convergence [68, 69], and 
it was concluded that solutions of bending problems of curved structures 
employing these elements do not always converge. In addition, since the 
use of flat elements for the solution of problems involving curved structures 
requires very large matrices, flat elements should be limited to problems 
involving arbitrary shells or shell cutouts, for which large matrices are 
obtained no matter what element is used. 

The inadequacies of flat elements led to the development of curved 
quadrilateral and triangular elements. The more general curved 
triangular elements can readily describe arbitrary cutout boundaries. 

The earliest constant curvature quadrilateral element was developed 
by Gallagher [70] by applying the shell theory (non-shallow) presented by 
Novozhilov [71], Although the displacement function chosen does not include 
rigid body terms nor does it lead to displacement continuity at the element 
boundaries, the results obtained converge satisfactorily. Geometrically 
more flexible quadrilateral elements were developed using shallow shell 
theory [72] and finite differences by Szilard and West [73], and subsequently 
by Tsui, et. al. [74], who also included the effects of shear deformation. 
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Bogner, et, al. [75] developed a cylindrical element wherein the dis- 
placement function includes both rigid body motion and satisfies the 
compatibility requirements. This element was later used to study non- 
linear shell behavior [76], More recently, a large variety of curved 
quadrilaterals have been developed: Conner and Brebbia [77, 78] (the 
Marquerre [79] shallow shell theory with nonlinear effects was employed 
in Reference 78), Cantin and Clough [ 80] (cylindrical element with a 
thorough discussion of the requirements for inclusion of rigid body motion 
terms in the displacement function), Wempner, et. al. [81] (transverse 
shear deformation is included), Ahmad, et. al. [82] (transverse shear 
deformation is included), and Key and Beisinger [ 83] (the displacement 
function is represented by Hermitian polynomials, and shear deformation 
is included). 1 

In Reference 73, vibration problems were first solved using curved 
quadrilaterals. More recently, Olson and Lindberg [84] established mode 
shapes and frequencies of curved fan blades employing curved quadrilateral 
elements. Greene, et. al. [8 5] also used a quadrilateral to study shell 
vibration problems. The shell theory utilized is the non-shallow shell theory 
of Novozhilov. The element of Reference 70 was used by Gallagher and 
Yang [86] in problems of stability of shells. 

Early work in the development of curved triangular elements was 
performed by Utku [87, 88], Prince [89], and Svalbonas [90]. Utku 
utilized the Marguerre shallow shell theory, and included shear deformation 
effects. Prince employed the three sub-element technique [91] to develop 
a constant curvature trianglar shell element, based on the non-shallow shell 
theory of Novozhilov. The three sub-element method Was also utilized in 
Reference 90 in conjunction with the non-shallow shell theory of 
Novozhilov, to develop a family of orthotropic, arbitrarily curved, triangular 
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shell elements. Curved triangular elements were recently developed by: 

Argyris and Scharf [92] (63 degree of freedom element in a ’'natural” 
coordinate system), Dhatt [93] (shallow shell theory with shear deformation), 
Strickland and Loden [94] (Novozhilov shallow shell theory), and Bonnes, 
et . al. [95] (three sub-domain method with Reissner [ 96] shallow shell 
theory). 

It should be noted, that the accuracy of solutions obtained by using discrete 
doubly-curved finite elements has mainly been verified only for simple classical 
problems. Thus, before definite conclusions can be drawn as to the 
suitability of some of these elements for solving problems involving shells with 

complex cutouts, additional, broader, comparisons must be made. 

Numerical Integration Method: This method was utilized in a general form 
in References 97, 98, and 99. Cohen [97] analyzed orthotropic shells of 
revolution, using the non- shallow shell theory of Novozhilov and the Runge- 
Kutta method of forward integration. The resulting system of equations is 
solved by Gaussian elimination. Both thermal and mechanical loads are 
considered, however, shell branching is not included. The unsymmetric 
loadings are expanded in Fourier series, and each harmonic is analyzed 
separately. Kalnins [ 98] solved similar isotropic shell problems using 
the Reissner shell theory, the Adams integration method, and the Gaussian 
elimination technique. Again the analysis is limited to problems involving 
two end boundaries. Mason et. al. [99] analyzed isotropic shells with 
completely arbitrary branching characteristics, subjected to axisymmetric 
load, using the nonlinear Love -Reissner -Kempner [ 100] shell theory and 
the Runge-Kutta method of forward integration. The arbitrary branching 
was accomplished with a finite element type (direct stiffness method) solution 
of the matrix equations. The nonlinear solution was identical to that of 
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Reference 10, discussed previously. Rung, et. al. [101,102] extended 
this method to the analysis of shells of revolution subjected to unsymmetric 
loading, using a linear shell theory. The direct stiffness matrix solution 
technique was adapted by Svalbonas [103] to orthotropic, second-order shell 
theory, including the effects of shear deformation and thickness stretch. 

In this reference, the nonlinear analysis of shells of revolution under 
unsymmetric loading is discussed, using a similar technique to that later 
utilized by Ball [ 20] (see finite difference discussion). The aforementioned 
method was applied to the analysis of shells of revolution comprised of 
orthotropic layers [104], and of sandwich construction with shear deformable 
cores [105], Another nonlinear analysis of shells of revolution, with two 
end boundaries, under an axisymmetric load, was presented by Kalnins and 
Lestingi [ 106] . In this reference, a form of Newton's method was employed 
to solve the nonlinear problem. The numerical integration technique, in 
connection with Guyan [152] reduction, was applied to the analysis of orthotropic 
stiffened shells subjected, to static loads in the STABS -II program developed for 
NASA by Svalbonas, et al. £107-110]]. This program may be employed in the solution 
of problems of shells of revolution with arbitrary geometry as shown in Fig. 1. 
Moreover, this program may be used in conjunction with discrete finite element 
analysis [m] . 

Calculations of natural frequencies and modes of vibration of shells of 
revolution using numerical integration was first accomplished by Kalnins 
[112], who included the effect of rotatory inertia. The problems solved 
entail isotropic shells wi+h two end boundaries. The frequencies are 
established by evaluating the frequency determinant. Orthotropic, ring -stiffened 
shells were analyzed by Cohen [113], who used numerical integration and a Stodola 
type [114] iteration technique. This iteration technique commences by assuming a 
value for the displacement components, setting the frequency to unity, and evaluating 
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numerically the inertia terms. These values of the inertia terms are 
substituted in the equations of motion resulting in a set of nonhomogene ous 
equations, the solution of which yields a first estimate of the mode shapes. 

An estimate of the frequency is then obtained by evaluating the Rayleigh 
quotient. This value of the frequency, together with the estimated values 
of the mode shapes are used to obtain a new set of numerical values for the 
inertia terms. The process then continues until the mode shapes obtained 
from two successive iterations vary by an acceptable error. It can be 
proven that this method, also known as the inverse power method [115], 
converges to the smallest frequency. 

The frequency determinant evaluation method as well as the Stodola 
method have certain advantages and disadvantages [116].. In the Stodola 
method, the lowest eigenvalue cannot be skipped, as is possible with the 
determinant evaluation method; however, if higher values of frequency are 
needed, the Stodola method requires modifications for "sweeping-out" the 
lower frequencies, and eigenvalue shifts to avoid slow convergence [ 114]. 

A generalization of the Stodola method was used by Cohen [ 117] to 
analyze the stability of orthotropic, ring -stiffened shells of revolution with 
two end boundaries, subjected to axisymmetric loading, employing the non- 
shallow shell theory of Novozhilov. Numerical integration by the Runge-Kutta 
method, and Gaussian elimination are used in the solution. The prebuckled 
state is obtained on the basis of a nonlinear solution and thus, the critical 
load and mode at buckling are established by solving a sequence of modified 
eigenvalue problems. The lowest buckling load corresponding to each 
Fourier harmonic buckling mode is established. The critical load at buckling 
is the smallest of these loads. A similar analysis including the capability 
of analyzing shells with "Y" branches, using the determinant evaluation 
method, was presented by Kalnins [ 118] . In this reference an attempt 



is made at establishing the critical load at buckling for shells of revolution 
of certain geometries subjected to single -harmonic un symmetric loading. 

Comparison of Methods ; The three methods discussed herein are approximate 
methods and, consequently, must be checked for accuracy or convergence 
of their results. The finite element and finite difference methods require 
at least two analyses with different grids to establish whether of the size 
of the grid used yields satisfactory solutions. The results obtained from a 
single solution by the numerical integration method, however, may be checked 
automatically for each shell segment. Moreover, the representation of a 
shell by finite elements involves geometric approximations which are not 
required in the numerical integration method. A disadvantage of the 
numerical integration method is that accuracy is lost when the shell is long. 

This however, is overcome by segmenting the shell into shorter pieces. 

A serious difficulty with the finite difference method is the instability 
of the solution at fine mesh sizes, and the slow convergence of the results 
obtained from single precision computer programs. For a stiffened cylinder, 
an example of the instability of the solution is shown in Figure 2. Introduction 
of double precision requires a reduction of the number of mesh points which 
appreciably limits the scope of a program [31] . It should be noted that the solution 
may not always converge asymptotically as the number of finite difference 
stations increases [40], In References 26, 27, and 119 it was indicated 
that the error associated with finite difference approximations of a given 
order is larger close to the boundary. Thus, the overall accuracy of the 
results may be increased if finite difference approximations of higher order 
are used at the end boundaries. Generally, in order to accurately establish 
the stress distribution near the shell boundary or at regions of high 
stress variation, or to properly locate shell discontinuities in curvature, 
thickness, etc., a mesh of variable spacing is required. However, the regions 
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Figure 2 Critical Axial Load as Function of Mesh Size for Cylinder Subjected 
to Axial Compression [From Reference 31 ] 
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of high stress variation cannot be predicted apriori and, consequently, 
the mesh of the first trial may be extremely fine in some regions 
of the shell and very coarse in others. In the use of variable mesh, the 
variation of the order of the error in the formulas, must also be considered. 

In the numerical integration method, however, the aforementioned difficulties 
are eliminated by automatically controlling the integration interval so as to 
obtain a solution of uniform accuracy. Finally, the errors in, the Runge- 
Kutta integration; formulas are of the order of h r whereas, in the finite 

differences employed in the aforementioned references, the error is of 

2 " • 

the order of h . In Figure 3 a comparison is presented of the results obtained 
by finite differences and by the numerical integration method for spherical 
caps subjected to uniform pressure. 

Finite element analysis of shells involves both mathematical 
approximations (those associated with a Rayleigh -Ritz analysis) and geometric 
approximations. The geometric approximations associated with revolved 
finite elements are discussed in detail in Reference 48. As previously 
noted, an advantage of the finite element method is that it may be 
employed to analyze arbitrarily branched shells in a routine manner. 

Numerical integration methods can also be employed to analyze arbitrarily 
branched shells [ 107], moreover, to obtain the same accuracy in the solution, 
much coarser idealizations can be used with the numerical integration method than 
with the finite element method. Finally, the components of stress and dis- 
placement obtained by the numerical integration method are of the same 
accuracy, whereas, the order of accuracy of the components of stress 
obtained by the finite element method is less than the order of accuracy of 
the components of displacements. 
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CHAPTER 1 

FORMULATION OF THE GENERAL NONLINEAR SHELL EQUATIONS 

Strain- Displacement Relations: The nonlinear strain- displacement 

equations for the Love- Re i s sne r - Kempne r shell theory are developed in 
Reference 100. A synopsis of this development ensues. 


Consider a deformable body in a state of stress due to surface trac- 
tions and body forces. The deformation in the neighborhood of a material 
point is defined by the unit elongations E a , Eg, E^, and the unit shears r a0’ 
Igy, r referred to a system of orthogonal curvilinear coordinates a, 0, y. 

(i=a, 0, y) represent the change of length per unit length, due to the deforma- 

til 

tion of a line element which was in the i direction prior to deformation. 

r. .(i, j = a, 0, y) represent the change in angle, due to the deformation between 
U ^ ^ 

two line elements which prior to deformation were in the i and in the j 

directions (i ^j). The unit elongations and shears are related to the strain 

components e..(i, j = a, 0, y) by the following relations [121] 
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The strain components, e_, are sufficient and convenient measures of the 
deformation in the neighborhood of a material point, defined in terms of the 
displacement components by the following relations [121] 
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The quantities ,©. . and oo . (i, j = a, p, y) are defined by [ 121] 
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u(ar, p, Y)» v(ar, p, ^f) and w(flf, p, y) are the components of displacement along 
the coordinates a, ft, y> respectively; w .(<*, p, y ) are referred to as the 
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components of rotation; H , HL , H , are the Lame coefficients of the 

s a . p V 

a, p, Y coordinate system. We shall make the assumption that the unit 
elongations, the unit shears and the rotations are small as compared to 
unity. On this basis, it can be shown [ 121 ] that 
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Taking into account that uf may be of the order of magnitude of <?-» 
Equations ( 1 - 2 ) reduce to 
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( 1 - 5 ) 


If the assumption is made that w. are of the same order of magnitude as 


, the above relations reduce to 


= e.. 
ij 13 


i» j = a, p, Y 


( 1 - 6 ) 


Consider thin shell. The position of points on the reference sur- 
face of the shell will be determined by the curvilinear coordinates a and 
P which are lines of principal curvature of the reference surface. The 
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position of a general point P of the shell will be specified by the coordinates 
a and p of the base of the perpendicular from point P to the reference 
surface, and the distance £ measured along this perpendicular. Referring 
to Figure 4, for a shell of revolution, these coordinates are designated by 
0, (p, and t . The Lame coefficients for a general shell of revolution may 
be written as 

H tt ar o< l4/r 2 ) 

H * r 4 (l-C/rj) (1-7) 

H = Hy = 1 
V t 

Referring to Figure 4, the Lame coefficients for the reference surface, when 
C/r « 1, are the radii of curvature (r Q , r^) of this surface, and must 
satisfy the Gauss -Codazzi compatibility relations [121] given by 

r o»<P = r i 003 * (1 " 8) 

In obtaining the above equations the following geometrical relation has been 
employed. 

r 0 = r 2 8 * n (*' ~9) 


In general, the displacement components in a shell may be expanded 
in a power series expansion of the t coordinate 


u(or> p,C ) = u(a, p, 0) + ^ + 
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Figure 4 Shell Geometry and Displacements 
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We shall assume that it is sufficient to retain only the first two terms in 
the above series. This is equivalent to the Kirchoff assumption that a 
line normal to the reference surface remains straight subsequent to de- 
formation. Moreover, we shall make the second Kirchoff assumption that r 
a line normal to the reference surface remains normal to the deformed re- * 

= e r - 0 at t - 0. In addition, we shall 
art 

assume that the second term in the expansion for w(«, |3, £) is small as 
compared to the first term and it can be disregarded. On this basis; Equa- 
tions (i -10) reduce for a shell of revolution to 

u(0, (p,i ) = u(0, (p) + 

v(0,<p,C) = v(0, <p) - Cw q (1-11) 

w(0,<p,C ) = w (0,<p) 

i r 

Substituting Equations (1-11) into Equations (1-3), and using Equations 
(1 -8) and (1 -9) it can be shown that 
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The above relations correspond to the Flugge-Byrne shell theory [12?., 123], 

They may be further simplified by disregarding — (i = 1, 2) as compared to 

r i 

unity. Thus, Equations (1-12) reduce to the following relations of 
Love-Reissner-Kempner accuracy c 
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(1-14) 


and Wq> w have been defined in Equations (1-12). 

For thin -walled shells, generally the rotation component is 

considerably smaller than the rotation components Wg, and may be 
disregarded. Thus, for thin shells, the non-linear strain Equations (1-5) 
can be further simplified to give 
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e ^Y = S cpC an< ^ ®ay = e 0£ are assume< ^ negligible compared to Sq^, 

Stress Resultant- Strain Relations : The stress resultants may be defined as 
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where T is the temperature distribution in the C direction, and the integrals 
are taken over the thickness of the shell (see Figures 5 and 6). In these 
definitions, the coordinate "C" is assumed negligible compared to "r". To 
this order of approximation, the difference between the inside and outside 
dimensions of an element of the shell becomes negligible [124]. If these 
assumptions were not made, then Nq^ and Mq^ 4 -M^q, since the radii 
associated with the cp and 9 directions are generally, unequal (r ^ 4 re- 
introducing the stress-strain relations for an orthotropic body in a 

state of plane stress into Equations (1-16), (1-17), using Equations (1-13), and 
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assuming that the reference surface is the centroidal surface of the shell, 
we obtain: 
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(K^j. Z> 3 ^ ) stiffnesses are defined as 
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Additional relations for different wall cross -sections are presented in 
Figure 7. 

Equilibrium Equations: In Reference 100 the following nonlinear 
stress equilibrium equations are obtained 
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Figure 7. Shell Section Properties 
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where F.(i = 0, cp, C) are the applied forces tangential and normal to the de- 
formed shell surface, whereas f^(i=0,cp t C) are the forces along the unde- 
formed coordinate system (see Ref. 137). The first three of Equations (1-20) 
are obtained by setting to zero the sum of the 9 , cp and C components of all 
the forces acting on a shell element. The last three of Equations (1-20) 
are obtained by setting to zero the sum of the moments about theoQ.i^, and 
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axis. As a result of the as sumption that C « r, the sixth equilibrium equation can- 
not be satisfied except for the special case of a sphere. However, within the frame- 
work of the present theory, this equation will not be employed in the solution of shell 
problems. 

Boundary Conditions: As shown in Reference 100 , for a unique 
solution either displacements or corresponding stress resultants may be 
specified on the boundary <p = constant. 
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The quantities T^q , 
and are defined by: 
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J , are referred to as the effective stress resultants 
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Convenient Form of the Shell Equations : In this section, the stress - 
resultant, strain relations and the stress -resultant equilibrium equations will be 
combined to obtain differential equations suitable for the Runge-Kutta integration 
procedure to be used in their solution [125-1281 • Be eliminating the strains 
and curvatures in Equations (1-18), using Equations (1-13) through (1-15) and 
the relation between the elastic constants v ^ Eg = v g. E^ , the following 
relations between the stress resultants and the displacements may be obtained, 
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Using Equations (1 -12£), (l-24a), and (1-2.5.C) 
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The final form of the differential equations necessary for the RUnge -Kutta 

integration procedure [128] may be divided into two groups. The first 

group of four differential equations is obtained by eliminating Q 0 from the 
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equilibrium equations, and by introducing the effective stress resultants 
defined by Equations (1-24) and their derivatives with respect to 0 , 
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The second group of four differential equations is" obtained by combining 
Equations (l-12e), (l-24a) and (l-25b, c, e). 
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In Equations (1-27) and (1-28) the stress resultants Ng, N^g,the resultant 
moments Mg, M^g , the rotation, , and the effective stress resultant, 
J , may be eliminated by using Equations (l-25a, d), (1-26), (l-24a, c) and 
(l-12f). The equations for these variables may be rewritten as, 
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Equations (1-27) through (1-29) represent a complete formulation of 
the nonlinear problem for a thin orthotropic shell, on the basis of the 
Love-Reissner -Kempner theory. Analogous formulations may be obtained 
by employing other shell theories. The basic differences between the 
various formulations will be in the coefficients of the equations, and in the 
number of differential equations. For example, in theories involving shear 
deformation ten basic differential equations may be obtained [ 105], 

For general doubly -curved Shells, the above equations may be 
written with the angle <p as the independent variable, whereas, for cones 
and cylinders they may be written more conveniently with the arc length 
s (s = rjdij p) as the independent variable. ; The stress resultants and dis- 
placements involved in these equations may be expanded in Fourier series 
in the 0 direction. Thus, Equations (1-27) and (1-28) will constitute a 
basic system of 8 first-order ordinary differential equations in the variable 
<p , and Equations (1-29) will constitute 6 algebraic equations. Notice, 
that derivatives of the shell geometric parameters do not appear in the coef- 
ficients of these equations. Moreover, notice that the Sjunknown variables 
are the quantities which enter into the appropriate boundary, conditions on 
the edge W = constant of a shell of .revolution. .... , .. .. • .. ... , 

Equations (1-27) through (1-29) will be solved by a forward numerical 
integration procedure in conjunction with the direct stiffness rqatrix mfethod 
(see Chapter 3). , The stress-. resultants and Q 0 , not involve d in the 

above formulation, may be obtained from Equations (l-24b) and (1 -20e), 
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•which may be rewritten as 
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The simplicity of the aforegoing formulation of the basic shell equations 
results in greater accuracy in the numerical solution. Notice that Equation 
(1 -30b) involves derivatives of shell geometric parameters. However, the 


computation of Q 0 is 4 secondary operation in the solution of the problem. 

In solutions of nonlinear unsymmetric loading problems with finite differences, 
this formulation has been found to yield more satisfactory results than one 
involving a basic system of four second -order differential equations [ 22 ]. 
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As previously mentioned, Equations (1-27) through (1-29) constitute 
a complete formulation of the problem for a homogeneous orthotropic shell. 
For eccentrically reinforced shells, this formulation must be revised. If 
we consider the reinforcement as being smeared over its spacing length, 

t 

a revision would be necessary in the Equations (1-18) to take into account 
the geometrical orthotropy and reinforcement eccentricity. This revision 
would affect only four of the Equations ( (l-28b, d) and (.1-29 a, b) ) if an 
appropriate shell reference surface is chosen. 

The revised integrated Hooke’s Laws ( Equations (1-18)) are derived 
for several cases of stiffening in Appendix A, Using these revisions, the 
following equations analogous to Equations (l-28b, d) and (l-29a, b) are 
obtained for shells with ring -stringer reinforcement: 
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•where the K. ., C.. and D. . stiffnesses are defined by Equations (A - 8 ) 
in Appendix A. Equations ( 1-31) have been derived on the basis of the 
stress resultant - strain relations (A - 7) in Appendix A. The following 
more general form of Equations (1-31), valid for multilayered shells with 
general ring and stringer reinforcement, may be obtained by employing the 
stress resultant-strain relations (A -9) in Appendix A. 


39 



^ - \ 4 “0 + ( K 22 + E^ Y 1 1 V N V ( V M V- (K 12 + -T§- 5 {?<"' 


. ,,12 , C 25 D 12 „ , >*• es*"' 8 sin, P , “9 

*1 V COS ip— VI SITHp) +JU - ( 5 Cjg) | 5 - + — cosy 

J ** r o ° 


“e, 0 _ 4- , K 22 D 22 

r l = ( 25 


K 22 C l5jl 


< V M T^-V N V (K i 2 - ~C^ _) F 0 (U ’ 0 +v cos( P 


1 2 

T ^ _ - 

2 <p* 

- 


-wsln 0 )-ljo -< C 15 +- C 


K 22 D 12, r w ’ 6S +U ’ 8 Sin,> "o “8 


+ COS (p 

o 


V<V*W ( k : 


K 12 C 25 r 1 L K 22 
K__ ' i 5 T 25 " C 


il-,. 

25 J 


N +( -? — 

T0 ( K 22 


- c i 5 )( c 


K 12 ,,, K 12 C 25 


K 22 D 22 s' 1 / 


K ?? D ??\- 1 K ?? K l ? K „ D . - l 

+ -^fj C7 5 < M 0 +M T^ < K n- - C 15>< C 25 + ^^f > < K 12 


-- Q 15 )|{~ (U^vcos^w sin</»+i 

25 J I o 


1 2lJ, K 12 C 15 

J r^' 


K 12 C 25 

C 14 M-tM - 5 ( 1 - 32 ) 


C 15)( C 25 + -^ 


25 y 


15 c ^lk 2 


fCl , + W^ + .^ + I 


e 

COS <0 

r ^ 

o 


2 \ -1 


M 9=<V M T0> -^ + (C 15 --4| 


K 22 + D„ 


£> +N T<0 )(C 15 


4o 




The above equations may replace the corresponding relations in Equations 
(1-28) and (1-29) which m^y then be combined with Equations (1-27) to form 
the complete set of equations for the analysis of problems involving a broad 
range of reinforced shells of revolution. In this formulation, the structure 
is symmetric about the axis of revolution, and thus, smearing for stringer 
reinforcement is Unavoidable unless the formulation is further complicated 
by expanding the circumferential stiffness in a Fourier series, as shown in 
Reference 9. Ring reinforcement properties, however, need not nec- 
essarily be smeared in such an analysis (see Appendix B), Indeed, in 
buckling problems, smearing of the ring or stringer reinforcing yields un- 
satisfactory results in cases where the half wavelength of the axial or the 
circumferential buckle pattern is smaller than the spacing of the reinforce- 
ment. The error resulting from smearing of the ring or stringer reinforcing 
for the case wherein the reinforcement of cylinders is spaced at exactly 
one -half wavelength of the buckle pattern, is approximated in Reference [ 120] , 
In order to eliminate this difficulty for ring reinforcement, discrete ring 
equations are obtained in Appendix B, and cast into a form suitable for 
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inclusion in the numerical procedure to be employed. These equations may then 
be utilized in cases in which the smearing technique is not suitable. 

It should be noted that Equations (1-27) through (1-29) do not apply 
at a closed apex of a shell. At the apex, the radius of revolution, r Q , 
vanishes, resulting in a singularity in these equations, as discussed in 
Reference 6 . This problem, however, may be circumvented as suggested 
in References 8 and 20 . The necessary differential equations and apex 
boundary conditions are derived in Appendix C. 



CHAPTER 2 


FOURIER SERIES EXPANSIONS 


Efficient techniques are not readily available for the numerical 

solution of partial differential equations of the complexity of those formulated 
in Chapter 1. However, by expanding the applied loading and the shell functions 

in Fourier series expansions in the cylindrical coordinate , 6, "it is possible 

to reduce the set of partial differential Equations ( 1-27, 28, 29) to sets of 

ordinary differential equations. The actual number of sets of ordinary 

differential equations will depend upon the type of load distribution considered, 

and the degree of accuracy required. These sets of ordinary differential equations 

may be solved by employing a standard numerical integration procedure such as. the 

Runge-Kutta [ 125 - 128 }. 

It is assumed that the applied surface loads can be satisfactorily 
represented by the terms for n = 0 , 1 , . . . N of their Fourier series expansion 
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Moreover, it is assumed that the displacement components, the rotations, 
and the stress resultants can also be represented satisfactorily with the 
terms for n = 0, 1, . . . N of their Fourier series expansion. 
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Linear Stress Analysis ; In problems involving linear stress analysis, 
or stability or vibrations of shells subjected to axi symmetric prestress loads, 
the substitution of the Fourier series expansions (2-1, 2) into the sets of 
partial differential equations (1-27, 28, 29) results in uncoupled sets of 
ordinary differential equations. These sets may be solved separately in 
establishing the amplitudes of the Fourier series expansions, which may 
then be employed in Equations (2-2) to yield the stress resultants and 
displacements. For instance, in problems of linear stress analysis of 
homogeneous orthotropic shells, in which the applied surface loads can 
be satisfactorily represented solely by the terms of the Fourier series 
expansion (2-1) having primed amplitudes, Equations (1-27) through (1-29), 
yield (N+l) sets, of the following relations — one set for each value of n(n»0, 1. . N) 
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N is determined by the accuracy requirements of the load representation 
(2-1) and of the solution. A similar set of equations may be obtained 
for surface loads described satisfactorily by the terms of their Fourier 
series expansion (2-1) having double primed amplitudes, by substituting 
the second portion of Equations (2-2), into (1-27) through (1-29), or by 
substituting -n for n in Equations (2-3). 

If a reinforced or a laminated shell is to be analyzed, the Fourier 
series expansions (2-1, 2) must be substituted into Equations (1-31) or (1 - 
instead of into the corresponding Equations (1-28) through (1-29). 

It should be noted, that in all the above mentioned cases, the axi- 
symmetric torsional case (n= 0) is uncoupled from the axisymmetric non- 
tor sional case. 



Nonlinear Stress Analysis: Nonlinear stress analysis problems for 
shells may be classified in two major categories, characterized by 
axisymmetric and by unsymmetric loadings. The problem of stress analys 
of an orthotropic homogeneous shell of revolution subjected to axisymmetri 
loading is described by the following equations obtained by substituting the 
Fourier series expansions (2-1, 2) with n= 0 into the set of partial differenti 


equations (1-27) through (1-29). 
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The nature of the nonlinear problem is evident from Equations (2-4) and 
(2-5). For example in a linear analysis, if solutions were obtained for a 
load described by Fq^ , and for a load described by F^^ , the sum of 
these solutions wmld represent the solution under a load described by 
(Fq + F^^ ). In a nonlinear analysis, the sum of the two solutions will 

not represent the solution for a load described by (F Q ^°^ + F r ^). 

0 b 

The presence of the (p derivatives in the right-hand side of Equations 
(2-5) does not result in additional complications. These derivatives could 
be eliminated by using Equations (2-4e, f). 

If reinforced or laminated shells are to be analyzed, the Fourier 
series expansion (2-1, 2), with n= 0, must be substituted into Equations (1-31) 
or (1-32) instead of into the corresponding Equations (1-28) through (1-29). 

The formulation and solution of the nonlinear problem under unsymmetric 
loading is more complex. In all prior formulations of this problem (20,22, 
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103, 104] , the applied loading could be described by a Fourier half- 
series. Thus, a line of symmetry was assumed in the loading distribution. 
If this assumption is not made, the full series must be employed. In this 
case, substitution of Equations (2-1) and (2-2) into Equations (1-27) through 
(1-29) will yield linear terms of the following form 

( f a' (a| CO. » 8 + f A^'inne) (2-6.) 

\n= 0 n= 1 / . / 

and nonlinear terms of the following form 
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where A and B denote the amplitudes of the Fourier series expansion. In 
order to eliminate the coordinate 9 from the equations containing terms of 
the form given by (2. 6b), the double series product terms must be con- 
verted to the form (2- 6a). This may be accomplished by using trigonometric 
angle difference formulae 
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Similarly, the last term of the right hand side of Equation (2. 6b) can be con- 
verted to a cosine series, whereas, the other terms of Equation (2. 6b) can 
be converted to a sine series. 

The nonlinear problem involving a homogeneous, orthotropic shell 
of revolution under unsymmetric loading can be described by substituting 
the complete Fourier series expansion (2-1, 2) into Equations (1-27) to (1-29). 
Then each of those could be reduced to the form 

[coefficient #1] cos n8 + [coefficient #2] sin n0 = 0 (2- 7b) 

The requirement that the coefficients of the cos n9 and sin n9 terms vanish 
simultaneously, yields two sets of equations. The first set is: 
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A similar second set of equations can be obtained from the vanishing of the 
coefficients of the sin n0 terms in the aforementioned equations of the form 
(2-7b). This set of equations involves the ( M ) harmonics. It can also be 
obtained by setting n to -n in the set of Equations (2-8), Notice that the 
harmonic amplitudes in each set are coupled. Moreover as can be seen from 
Equations (2-9a, b) both sets of equations are coupled through the non- 
linear terms. Thus Equations (2-8) cannot be solved separately for each 
value of n, as in the linear case. The nonlinear terms of Equation (2-8) are: 
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second set of equations (obtained by a setting -n for n and the double primes 
for primes in Equations (2-8) ) are given by 
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•where a= 0 for n= 0, and a = 1 for n ^.1. It is evident, that by assuming that a 
line of symmetry exists in the loading (as done in earlier references), only 
the symmetric or the antisymmetric half of the Fourier series expansion (2-1, 2) 
can be used, and consequently only one of the two aforementioned sets of 
equations is required. If we choose to use Equations (2-8) and (2-9a), for 

example, the latter are simplified further since all the double primed terms 
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become zero. As noted previously, revised equations can be formulated 
for the reinforced or laminated cases. 

As will be discussed in Chapter 3, Equations (2-3) or (2-4) and (2-8) 

will be employed to establish the stiffness matrices and the load matrices 

of the various segments of the shell. The numerical solution of the equations 

for the linear analysis (2-3) may be attained by a technique described in 

Chapter 3. To apply this technique to the non-linear equations they must first 

be linearized by the use of a suitable method. For instance in Reference 20, 

the harmonic amplitudes of the shell functions are established throughout the 

shell by first solving the unsymmetric linear problem. The results may be 

employed in Equations (2-9a, b) to establish a numerical value for the nonlinear 

terms. In general, the linear analysis may yield a number of zero harmonic 

amplitudes. For instance, the linear analysis will yield nonzero harmonic 

amplitudes only for the values of n corresponding to a non- zero load in 

Equations (2-1). Thus, if for example, the load is described by the n= 1, 2 

( 1 ) ( 2 ) 

harmonics, the linear analysis will yield only the Qq ' and Qg ' harmonics 
of (JUq. However, to establish the values of the non-linear terms in Equations 
(2-9a, b) additional harmonics of OUg are required, as for instance fig"^ and 
fig 1 ); where n= 2, 3. . . P, P being the harmonic at which the infinite series 
is truncated. To obtain the values of the non-linear terms, all the harmonics 
of the functions which are not established by the preceding cycle of the 
analysis, are set to zero. The numerical values of the non- zero non-linear 
terms are computed, and subsequently employed in Equations (2-4) or (2-8) 
resulting in an uncoupled system of equations which may be solved to yield 
another set of values for the harmonic amplitudes. It should be noted, that 
some of the equations for the harmonic amplitudes (2-4) or (2-8) may not 
contain actual load terms, but rather the numerically represented non-linear 


.65 



terms. Subsequently, the values of the harmonic amplitudes are used to 
compute a new set of non-linear terms, and the process continues until the 
harmonic amplitudes are established to the desired accuracy. Moreover, 
it should be noted, that depending upon the number of non- zero nonlinear 
coupling terms established in each cycle of the analysis, the number of 
equations to be solved may vary per cycle. Although the aforementioned 
technique is the most straightforward, it may “not converge, however, for 
certain cases [22, 131], 

In another method [9, 14, 15, 99] the harmonic amplitudes are also 
established from the solution of the linear part of Equations (2-4) or (2-8). 
Instead of utilizing the numerical value of both terms in the non-linear products 
of Equations (2-9a, b) only one term is substituted numerically. Substitution 
of the linearized Equations (2-9a, b) into Equations (2-4) or (2-8), yields 
sets of coupled linear equations, which may be solved by any method that 
does not result in extremely large matrices; as for example, by a numerical 
integration method analogous to that presented in Reference 107. The fore- 
going methods are adequate only for problems in which the effect of non-linear 
terms is small [14, 131], although in Reference 20, it is indicated that these 
methods may be employed in the solution of problems involving relatively 
large deflections. 

In the solution of non-linear problems the load should be applied in 
increments* When the solution converges for one load level, the load is 
increased by another increment and the process continues until the components 
of stress and displacement of the shell are established for any value of the load 
desired. If the load is continually increased, it will reach a value for which 
the solution diverges. This indicates > that for this value of the load, the corre- 
sponding deformation pattern of the shell has reached a level of instability. 
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For ideally "perfect" shells having identical prebuckling and buckling defor- 
mation patterns, the classical bpckling loa^jis ^e value pf ^the load prior 
to that for which the solution diverges. ,. IJeiypypr, ideally ''jperfect'.' shells, 
whose buckling deformation pattern, differs. from their prebpckling deforma- 
tion pattern may atain one or more states of instability for smaller values 
of the loading, which will not be obtained by the analysis outlined above* 

Other methods for establishing all the states of instability are presented 
in Chapter 4. 

It should be indicated, that the af o r erne ntione d s olution of the non- 
linear equations may not converge for values of the load less than the maxi- 
mum value for which the prebuCkling deformation, pattern of the shell is 
stable [22, 41,131]. Thus, the aforementioned solution may yield conserva- 
tive values of the load for which the prebuckling deformation pattern becomes 
unstable. .. 


A method of proven convergence for the solution of nonlinear equa- 
tions is Newton's method [17,18], wherein each harmonic amplitude is ex- 
pressed as !the sum of two parts , ran assumed solution, and a correction to 
the assumed solutions- . >■ ’ !:' > • . "■■■) > 


Y . = Y + AY 
m + 1 m 


(2-10) 


Equation (2-10) may be substituted into either Equations (2-4, 5) or Equations 
(2-8, 9a). The resulting equations will contain terms of the type (Y m )(AY) 
as well as (AY)^. In these equations the Y^s are numerically known, from 
a previous iteration. In a buckling analysis the nonlinear (AY) terms may 

I . . . • v . > , ■ f . , .... , , "... ...... 

be neglected. For a postbuckling analysis, however, these terms should be 

retained [1’29, 130], Restricting ourselves to a nonlinear analysis for the 

prebuckling state, the linear AY correction equations will be uncoupled or 
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coupled depending upon whether axisymmetric (Equations (2-4) and (2-5)), 

or unsymmetric (Equations 2-8) and (2-9)) loading is considered. The 

procedure is as follows: First, the problem is solved for a small value 

of the load, where the linear theory is accurate. This may be accomplished 

by using the aforementioned AY equations, but setting all the terms to 

zero. The solutions for AY, and = 0, are substituted into Equation (2-10). 

The values of Y , thus obtained are Substituted for the Y (which in the 
nit i m 

previous step were set to zero) in the aforementioned AY equations, pro- 
ducing a set of linear, coupled (if the load is unsymmetric) equations. These 
equations are now solved for new AYs, and these, with the current Y^s, 
are substituted in Equations (2,-10) to produce a new set of Y m + ^s. This 
procedure continues until the harmonic amplitudes are established to the 
desired accuracy. The load on the shell is now increased by an increment 
and the whole procedure is repeated. The solution obtained for the previous 
value of the load, or a set of harmonic amplitudes obtained by extrapolation 
can be used as a starting point. Non convergence of the solution indicates 
that the last load increment has increased the value of the total load above 
the limit point (the point of zero slope of the load- displacement curve). The 
solution may be repeated using a smaller load increment, and values of the 

• • i 

total load as close to the limit point as desired may be obtained. The last 
value of the load for which the solution converges is taken as the maximum 
value of the load for which the prebuckling deformation pattern of the shell 
is stable. 

Reference 22 essentially utilizes the foregoing technique with one 
important simplification. The nonlinear terms which induce coupling of the 
harmonics are disregarded. This is necessary inasmuch as the finite dif- 
ference technique employed in Reference 22 results in large matrix equations. Con 

sequently, their solution becomes extremely time-consuming if the harmonic 
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are coupled. The penalty for this simplification is a slower convergence, 
and for more complex problems, possible divergence at loads less than the 
• maximum value of the load for which the prebqckling deformation pattern 
of the shell is stable. In a numerical integration procedure of the tyoe 
presented in Reference 107, the resulting matrices can be maintained at a 
reasonable size, and consequently it is no t t essential to make s implif cation s . 

Thus, by employing any of the aforementioned methods, the non- 
linear problem simplifies to that of successively solving linear equations, 
which in the most complex case (unsymmetric loads) involves harmonic 
coupling. The procedure for the solution of these types of equations Is 
presented in Chapter 3. 

■ 1 . -. ; ' ■' ■ : , ■ " 

For a postbuckling analysis, as noted earlier, the nonlinear terms in 

AY must be retained in order to cross bifurcation points [129, 130]. In addi- 

. ... I ... . ■ : • • 5 

tion, a method of proven convergence, such as the Newton iteration method 

must be used in the solution. Some difficulties do exist, however, in utilizing 

. . ■ • ■ ■ . ' t" ■ . ■ .1 

the equations for postbuckling investigations. The most important difficulty 
is exemplified by the case of a shell under axisymmetric loading only, wherein 
a postbuckled configuration may exist which is described by harmonic ampli- 
tudes other than the zeroth. In this case, a small "disturbance load" 
involving several harmonic amplitudes must be added to allow the shell to 
deform into the proper configuration in the postbuckling range. This dis- 
turbance load must be small in magnitude, as compared to the primary 
axisymmetric load, but it must also be described by a sufficient number of 
harmonics so that an adequate description of the postbuckled configuration 
may be obtained. 
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CHAPTER 3 


SOLUTION FOR STATIC RESPONSE 

Shell structures found in spacecraft, aircraft engines, or submersibles, 
are usually comprised of several, interconnected, singly and doubly curved 
shells* For a given shell geometry, the accuracy of the results obtained 
on the basis of the numerical integration method is dependent upon the length 
of the shell along the meridian. If the shell is long, the effect of the stress 

resultants and displacements at the one edge will have a negligible effect 
upon the stress resultants and displacements on the other edge. This could 
result in a number of terms in the stiffness matrices which are inaccurate, 
inasmuch as they constitute small differences of large numbers. This 
difficulty may be circumvented by subdividing the shell into segments by 
introducing fictitious boundaries. Such an approach is amenable to the use 
of local coordinate systems, and includes the unique self-checking features 
discussed in the Introduction. Equations for establishing the appropriate 
lengths of segments for shells of various geometries are presented in 
Reference 108, for linear and nonlinear analyses. 

In this chapter, we shall first obtain all the matrices pertaining to 
single shell segments.- Then we shall proceed to couple these matrices 
together and apply the boundary conditions, in order to obtain an overall 
matrix equation describing the equilibrum of the total structure. The 
solution of this equation yields the stress resultants and displacements, 
at the joints. These stress resultants and displacements are used in 

TO 





establishing the distribution of the stress resultants and displacements 
throughout all the shell segments. 

Coordinate Transformation Matrices: A. typical shell segment is presented 

in Figure 8, Various possible geometries of shell segments are given in 
Figures 9 through 13. Since the various shell segments may be described 
by different coordinate systems and different geometric variables, the 
stress resultants and displacements referred to the coordinates of a seg- 
ment must be transformed to the reference global coordinates. The edge 
forces on the typical shell segment in this global Z, R,0 system are shown 
in Figure 14. The components of stress resultants and moments referred 

i ■ * v 

to a local coordinate system are denoted by Greek subscripts(T^g, N^, M^), 

whereas the components of stress rsultants referred to the global coordinate 
system are denoted by Latin subscripts (F^, F^., F^). The global 
moment is denoted by M. The appropriate coordinate transformation 
matrices then are: 


{F(i)} = [IFT ] (f(i)} 
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Figure 12 Cylinder 
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8 


(F(j) } = [ JFT] (f(j)} 


and 


{ A (i)> =|IDT] {5(i)> 


(A(j)}=[JDT] {.6 (j)> 





where [ IFT] , [ JFT] , [ IDT] , [ JDT] denote the I Force Transformation, 
the J Force Transformation, the _I Displacement Transformation, and 
the Displacement Transformation matrices, respectively. The subscripts 
i and j refer to the meridional coordinates of the beginning and end of the 
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segment, respectively. These transformation relations are valid for the 


functions F(0,<p ) and £(9,<p ) as well as, for the amplitudes of harmonics, 

K \ r K 4 - 

F (n) (y k ) and £ (n) (<p k ) , (k = i,j). 

In the sequel, the pertinent matrix equations will be written for one 

. > ■■■ ■■■ ' ' ' V 

harmonic at a time; this will not result in any loss of generality of the 

t '.r 

solution to be discussed. However, the harmonics are coupled for non- 
linear problems having an unsymmetric load. Thus, it would not be possi- 
ble to write the matrix equations for only one harmonic. Consequently, the 
size of the matrices would be multiplied by N, the number of harmonics 
to be retained when the Fourier series expansions are truncated. Hence, 
for a nonlinear problem with unsymmetric load, if N harmonics are re- 
tained in the Fourier expansions, a typical transformation matrix may be 
denoted by [IFTN], and assumes the form 

[ IFT] 0 — 0 — 

0 \ 

I ^ 

0 [IFT] 

1 

0 ; . w» mm . 0 

where "N" qf the 4x4 [IFT] matrices are located on the diagonal. Notice, 
that for nonlinear problems involving unsymmetric loads, the other matrices 
such as the stiffness matrix, may not be block diagonal matrices. Such 
special matrices will be developed separately as the need arises. 

Segment Stiffness Matrices; The suitable differential equations for each 

specific shell segment are solved for different sets of initial conditions by 

the Runge-Kutta forward integration method. Any satisfactory Runge-Kutta 
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[IFTN] = 
4Nx4N 




formula may be used [ 125] - [ 127] . The one employed in this investigation 
is [128]; 


Y, = Y + Y 

1 m 2 m 


Y = Y + Y, 

2 m 2 1 


(3-6) 


Y = Y + At Y, 
3 m 2 


At * * 

Y = Y + ===- (Y + 2 Y, + 2 Y, + YJ 
m+1 m 6 m 1 2 3 


These relations may be employed to establish the value Y , of a 

m+ 1 

function at point (m+1) if the value Y of the function and its derivative 
Y ^ with respect to the integration variable, are known at point m . The 
symbols used in the above relations are defined as: 


At = the integration interval from m to (m+1). 

Y = the value of the function at point m obtained by 
numerical integration, 

Y^ + = the value of the function at point (m+ 1 ) obtained by 

numerical integration. 



Y^ = the first and second estimate respectively, of the 
values of the function at the mid-point between 
points m and (m+1) . 


Y = the first estimate of the value of the function at point 
3 (m+1) . 


O 

( 


) = derivative of the function with respect to the integration 
variable. 
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The Runge-Kutta integration method is employed in establishing the 

ifcll 

values of the functions at the j edge of a segment from the assumed values 

th. 

of the functions at the i edge of the segment. A number of points are 

automatically chosen along the meridian of the shell segment, 

i , i+1 , i+2 ... m . . . j-2 , j-1 , j . The spacing of the points is denoted 

by At , and may vary from point to point. The derivatives of the functions 
th 

at the i boundary are established initially by Equations (1-24) through 

til 

(1-26) . The values of the functions and their derivatives at the i bound- 


ary are then employed in establishing the value of the functions at point 

^ th 

(i+1) . The process is repeated until the values of the functions at the j 
edge are established 

The process for establishing the values of the functions ait point m+1 

from those at point m is as follows. The values of Y (stress resultants 

m 

and displacements) are employed in Equations (1-24) through (1-26) to es- 

• • 

tablish the derivatives Y . The values of Y and Y are employed 

m m m 

to compute Y^ . These are the values of the predicted stress resultants 
and displacements at the point midway between m and m+1. Subse- 
quently, they are employed in Equations (1-24) through (1-26) to establish 
the derivatives Y^ . The values of Y^ are then computed from Equation 
(3-6b) . The values of Y_ represent a corrected estimate of values of 
the stress resultants and displacements at the point midway between 
m and m+1 . These values are then used in Equations (1-24) through 
(1-26) to establish the derivatives Y^ * Subsequently, the values of Y^ 
are computed from Equation (3 -6c) and employed in computing the values of Y 
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The values represent the first estimate of the values of the predicted 

stress resultants and displacements at the point (m+1) . A corrected esti- 
mate ^ of these functions is then computed using Equation (3-6d). 

Subsequently, the values of Y^ and are compared, and if they 

differ by less than a set tolerance, d^ , the process is continued to es- 
tablish the values of the functions at a point (m+2) , located at an interval 
2 At from point (m) . If the values of Y„ and Y . differ by more 

than the set tolerance, d^ > the current At is halved, and the process is 

repeated until the values of the functions Y. and Y differ by less 

3 m+1 

than the tolerance dj . Using the same interval At, employed in the 
previous step, the process is continued to establish the values of the 
functions at point (m+2) . Thus, the interval At may vary from point to 
point. This procedure is referred to as automatic step control, and provides fo 


uniform accuracy in the solution of the differential equations. 

As indicated previously, we will start by assuming the values of the 

th 

eight shell functions (stress resultants and displacements) at the i 

boundary t>f each segment, and we shall compute the corresponding values 

of the eight shell functions at the boundary. To isolate the influence 

of each of the eight functions we Will set one function to unity and the 

others to zero. The process is outlined schematically in Fig. 15. In 

th 

columns 1 through 4, one of the displacements at the i edge is succes- 
sively set to unity, while the remaining displacements and the stress re- 


sultants are set to zero. The corresponding values of the harmonics 

/°^(j) of the stress resultants at the edge are represented by the 

__ (n) 

matrix X , whereas the values of the harmonics 6 (j) of the displacements 

8k 




Figure 15 Calculations for Influence Coefficient and Load Coefficient Matrices 
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are represented by the matrix V-.. In columns 5 through 8, one of the 

th 

stress resultants at the i edge is successively set to unity, while the 

remaining stress resultants and all the displacements are set to zero. 

The corresponding values of the harmonics f^(j) of the stress resultants 
"til 

at the j edge are represented by the matrix X~, whereas the values of the 

(n >. 9J 

harmonics 8 (j) of the displacements are represented by the matrix 

Ld 

th 

In columns 9 to 18, the stress resultants and displacements at the i edge 

(n) 

are all set to zero. The values of the harmonics f (j) of the stress re- 
th. 

sultants at the j edge due to the external distributed loads acting along 

the segment of the shell are represented by the matrix, whereas the 

/ n \ 

values of the harmonics S' (j) of the displacements are represented by the 
matrix V-y Notice, that as many loading conditions as desired can be 
considered. In Fig. 15, ten loading columns are shown. The matrices 
K , X, , X, , fl , u, % are referred to as the influence coefficient 

1 U $ 1 U D 

matrices. Each of the eight different edge conditions, and the ten load- 
ing conditions produce a column in the appropriate influence coefficient 


matrices. 


Fig. 15 is applicable to a single harmonic of a linear problem or to 
a nonlinear problem with an axisymmetric load [ 1 0 7l • The nonlinear 
problem associated with unsymmetric loading is more complex. Fig. 

16 represents a schematic diagram for establishing the influence coef- 
ficient matrices for nonlinear problems entailing, for example, the coupl- 
ing of harmonics n and n f . The initial matrices for the values of the func- 
tions at edge i (the initial conditions) and the influence coefficient matrices 
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are now 16 x 16 matrices. In Fig. 16, the star matrices denote the effects of 
harmonic coupling. If there were no coupling, these matrices would be null mafrric 

The differential equations for the nonlinear problem may be linearized 
by the Newtonian method. (As discussed in Chapter 2, the equations in 
A Y are linearized by dropping the nonlinear terms in AY , and by using 
the previous values of Y in the product terms (Y) ( AY)). Thus, super- 
position is possible for both the linear and the nonlinear problem. Con- 
sequently, the influence coefficients may be employed to yield the stress 

th 

resultants and the displacements at the j edge in terms of the actual stress 

til 

resultants and displacements at the i edge. Using Equations (3-1), to 
(3-4) , the stress resultants and displacements may be expressed directly 
in global coordinates as 


{ F ( j ) } = [JFT] (f(j)} = [JFT] pj :x z :X 3 ] 6 (i) 

i 


(3-7) 


{A(j)} = [JDT] { 6 ( j ) } = [JDT] [U x I \ 



(3-8) 


where {f } is a scaling factor for the load. If the loading cases employed 
in establishing % and ^ are the actual loadings considered, then the 
loading vector {j? } is a unit vector. Solving Equation (3-8) for vector 
{f (i) } , and employing Equations (3-3) and (3-1) to convert { 6 (i ) } and (f(i)} 
into global coordinates, we obtain 

(F(i)} = [IFT] [^ 2 ]" 1 A JDT 1 T UDT] T {A(i)} - {f } 
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Using Equations (3- 1) and (3- 3) to convert the vectors { 6 (i) } and (f(j)} in 
Equation (3-7) into global coordinates, and substituting in the resulting 
equation, the values of {F(i)} from Equations (3-9), we get 

(F(j)} = [JFT] ptj [lDT]{A(i)} + [JFT] [xj [g/J" */[ JDT] T (A(j)} - [^][IDT] T (A(i) } 
- [^ 3 } {f}j + [ JFT] pj {f} (3-10) 

Equations (3-9) and (3-10) may be combined in the form 


{F} = [k] {A} + {L} 

where, referring to Equations (3-9) and (3-10), we may write the stiffness 
and load matrices in a combined form: 



Equation (3-11) may be verified by carrying out the matrix multiplication 
and comparing the resulting matrix with Equations (3-9) and (3-10). As 
evident from Equation (3-11), in order to compute the stiffness matrix [k] 
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and the load matrix [L], it is necessary to invert only the matrix. 

(V-^ is a 4x4 matrix in uncoupled cases, 4N x 4N in the unsymmetrically 
loaded nonlinear case. ) 

In other numerical integration methods, [97, 98] the influence 
coefficient matrices X and $4 are used directly instead of calculating 
the stiffness matrices of the shell segment. Ther 6 are, however, many 
significant advantages in employing the direct stiffness technique. 

The procedure, employed in this investigation for the solution of 
the boundary value problem, subsequent to the formulation of the stiffness 
matrices of the shell segments, is exactly that employed in the finite ele- 
ment techniques. Thus, all the matrix manipulation methods developed 
for finite element solutions may be utilized in the present method, as for 
example, those for large scale matrix inversion, using packing techniques, 
or taking advantage of banding, etc. Moreover, the direct stiffness method 
can be applied without modification to arbitrarily branched shells, as well 
as to shells with discontinues changes in meridian slope. Finally, the 
other techniques [ 97, 98] are more efficiently utilized with the use of 
Gaussian elimination procedures, such as Potters' [ 4] procedure. These 
methods however, may be prone occasionally to error accumulation in the 
calculation sequence [ 132], whereas the Choleski method and the transfer 
matrix method utilized herein, tend to involve negligible error accumulatior 
Therefore, many finite difference schemes, utilizing some form of the 
Gaussian elimination technique require double precision arithmetic. 
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Properties of the Stiffness Matrix [ kl : In the solution of shell problems. 


for self- checking of the arithmetic , it is convenient to utilize the fully 




symmetric matrix [k] , defined by 


/x 

[kl* 


2ir r (i) 
o 


O o o © 


2ir r (j) 


Ik] 


(3-12) 


and therefore. 


/x 

[L] = 


N 2tr r (i) : 


X* 


! 2 ir r (j) 
° x 


[L] 


(3-13) 


so that 


/N /x 


[F] = [k] [A] + [Li] 


(3-14) 


where 


A 

[F] s 


F (i) 1 

1 

ss 

‘2ur o (i) F(i)" 

1 


i 

F(j) . 

j . 

_2irr o (j) F(j) . 


F is measured in units of force/unit length, and F is measured in units 
of force. 


In the case of the stiffness matrix [k], for either linear problems 
or nonlinear problems with axisymmetric loading, the matrix required to 
convert [k] into the symmetric matrix [k] may be obtained by inspection. 
This is not the case, however, for a nonlinear problem with unsymmetric 


91 



loading. There is no apriori reason in this case to assume that the stiff- 
ness matrix of the coupled harmonics can be converted to symmetric form 
inasmuch as this matrix relates a combination of harmonics of forces and 
displacements. Thus, the existence of a symmetric matrix must be proven. 
Consider, for instance, the case wherein the zeroth harmonic is coupled 
to any other harmonic, for example the In harmonic. For this case, the 
force-displacement relations may be written symbolically as : 


F (0) 

1 

- k (°,o) ; k (o, N> - 



f (N)' 

1 ’ 

k (N,0) : k (N,N)_ 

"■ j 

"Jm 


For the shells under consideration,, in the absence of body forces, the 
Betti -Rayleigh reciprocal relations [133] are (where the primed quantities 
belong to one system of forces and displacements, and the unprimed to 
another) : , 

i }*\ 

J F(0) A (6) d 9 = J F(0) A (0) d 0 „ (3-16) 

0 0 


Inasmuch as we assume that the zero harmonic is coupled only with the 
N** 1 harmonic, we may write 


(0) 

A(0) = A V + 


a (N) 

A cos 


N 0 


£( 0 ) 



+ A 


(N) 


cos N 0 


F (0) = F (0) + F (N) cos N0 

* «(0) '(N) 

F(0) = F F K ’cos N 0 


(3-17a) 


Substituting Equations (3- 17a) into Equations (3-16) we obtain 
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/( F (0) + F (N) cosN0) (A (0) + A (N) cos N9) d 0 = 

0 

2 IT ' 

f (f'(°K f'^cos N 0) (A (0 ^ + A (N) cos N 0 ) d 0 
0 

Integrating we obtain 

2„F«»a' ( 0) + n F (N) a' (N) = 2w f' < 0) A (0) + w f'< N >A (N) (3 ' 17b > 

The force -displacement relations for the two sets of harmonic amplitudes, 
on the basis of Equation (3-15) are given by 

F (N) = k (N>0) A (0) + k (N ’ N) ^ N) F , ( N ) = k (N,0)^'(0) + k (N,N) A '(N) 

F (0) = k (0,0) A (0) + k (0,N) ^(N) f '( 0) _ k (0,0) ^'(0) + fc (0,N) A '(N) 

Substituting the above relations into Equation (3- 17b) we get 
2 it k ( °’ 0) A <0) A (0) + 2* t («.V i '(0) t , 1; (K,Oyo) i 'TO tt k «W(»] 

= 2tt k< 0 * 0 >A (0, A (0) + 2ltk «,N) A '(N) A (0) + ^ k (N,0)^’(0)JN) +w k (N,N)^(N)JN) 
This relation may be rewritten as 

2 k (0 ' N) [A (N, A l(0) - a' (N) A ( 0) ] - k (N ’ 0) [A (N) A (0) - a' (N) A (0) ] 

This relation is valid if 

2 k<0’ N) _ k (N, 0) . (318) 

thus indicating that the stiffness matrix of Equation (3-15) is unsymmetric. 
However, it is relatively simple to form a matrix which can transform the 
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stiffness matrix of the coupled harmonics into symmetric form. Note, that 
each individual diagonal block (k^' ^ or k^’ ^) of the stiffness matrix 
for the coupled harmonics is not a symmetric matrix, but may be converted 
readily into a symmetric matrix (k^’ ^ or k^’ ^) by Equation (3-12) . Thus 
in this specific case, the appropriate symmetric matrix corresponding to Equation 
(3-15) is given by 




k (o, o) : k (o,N) 
k (N, 0)1 k (N/N) _ 


(3-19) 


The following notation is introduced in order to identify data in 
subsequent discussions and calculations: 


[ F] 


(n) _ 




k ] {n) U] 


(n) 




L] 


(n) 


(3-20) 


where: 

th 

s indicates the s segment of the shell connecting joints i and j, 

n indicates the Fourier harmonic. ( For coupled harmonic problems, 
the matrix superscript would be (n, n',n") where n, n, 1 n" are the 
coupled harmonics. ) 


Structure Matrices and Stiffness Analysis: The direct stiffness method 

[134] is employed in calculating the interaction of the segments comprising 

the shell structure. To increase the capacity of the program, the shell segments 




Figure 17 . Example of Region Topology 
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■will first be coupled into regions using a Guyan 0-52] reduction procedure. These 
regions are defined as singly-connected shells, with no internal concentrated line 
loadings (Figure 17). The next step is to construct the region stiffness matrix 
and the matrix of fixed-end forces W- This requires splitting each 
segments £k J matrix into its four 4x4 matrices (for coupled problems 4 Nx 4N ) , and 
inserting the portions into the region stiffness matrix in accordance with the 
topological arrangement. The [l] matrix is similarly split into two 4xP matrices, 
where P is the number of individual loadings considered separately. (For nonlinear 
cases, the stiffness matrix changes with the load, consequently, only one loading 
case can be considered at a time. Thus, the split load matrices are 4x1 for 
axisymmetric loading, and 4Nxl for the unsymmetric coupled problem. ) Thus, in 
addition to the geometric description of each segment, its position in the assembly 
must be specified. To this end, all segments begin (i) and end (j) at a joint. 

The s th segment is said to connect the i th and j^ h joints. (Not the j th and i th 
joints, since direction of increasing coordinate within the segment must be from i 
to j). To allow for the possibility of discontinuous centerlines within a region, 
kinematic links must be included. These links are rigid pieces which relate 
displacements across a discontinuity. Thus a kinematic link matrix £skl] must also 
be formed. Due to the topology and line-load requirements for regions, the equation 
of the coupled segments will be the following: 



( 3 - 21 ) 
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where 



SKL 



SKL n! SKl i 2 

I 

| 


ij^i : 


SKL 


22 



and where iR, jR refer to the region initial and final points, and the B M. 
and jjj'J are the deflection, stiffness, and load matrices of internal segments. 
If there are no internal kinematic links, £SKL J will be an identity matrix. 
Partitioning Equation (3-2l) will yield: 




N -WW-MH-m 

M -[£][■.] *[C]H*fl 


(3 -22a) 
(3 -22b) 


Solving Equation (3-22b) for [A J and substituting into Equation (3-22a) yields: 

ra - csm] + tQ] 

8xP 8x8 8xP 8xP (3-23) 

where 


/s / ,/\ S\ -1 /\ \ 

M = [ K ll] - ( [ K 12][ K 22] C K 2l] ) 

/V / /\ /\ /S, \ 

Chi] ■ [hJ - ( t K i2][ K 22 ] t L ] ) 

A 

The next step is to construct the total structure stiffness matrix j*K ^ "1 and the 

total structure load matrix ^ J . This requires the splitting of the stiffness 
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matrix J of each region into its four 4x4 component matrices (for coupled 


problems 4 Nx 4N ) , and inserting the portions into the total stiffness matrix, in 

•• • . - • v : : A ;; ■! 

accordance with the topological arrangement of the structure. The matrix is 
similarly split into two 4xP matrices, where P is the number of individual loadings 
considered separately. (For nonlinear cases, the stiffness matrix Changes with the 
load, consequently, only one loading case can be considered at a time. Thus the 
split load matrices are 4x1 for axisymmetric loading, and 4Wxl for the unsymmetric 
coupled problem). Therefore, in addition to the geometric description of each 
region, its position in the assembly must be specified. The initial point of all 
regions will be denoted by (i) whereas the end point will be denoted by ( j ) . 
Inasmuchas there are four degrees of freedom at each joint, for a shell with M 
joints the total stiffness matrix is 4 Mx4m ( 4MNx4MEf for coupled problems). Hence 
using equilibrium relations for all the joints, we can form the following 
equation 


[k] t [a]^ 


(3-24) 


The subscript T denotes a matrix which includes terms for all the joints. 
Equation (3-24) characterizes the structure without taking into account any 
external boundary conditions. For axisymmetric and antisymmetric (n=0, 1) 
load cases, the matrix [K],^, can be singular. This may be physically in- 
terpreted, as follows. The stiffness matrix permits calculation, of the 
stress resultants from the displacements; thus, the inverse of the stiff- 
ness matrix would relate displacements to the stress resultants. The 
displacements, however, are not unique inasmuch as one valid solution may 
differ from another by rigid body motion. Hence, it can not be anticipated 
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that a relationship may be established relating all the valid displacement 
fields to the unique set of stress resultants; indicating that the total stiff- 
ness matrix is not invertible, that is, it is singular. However, the total 
stiffness matrix of a complete shell of revolution for harmonics other than 
n*0, 1 need not be singular. For harmonics greater than unity, the stress 
resultant harmonics are self -equilibrating. Moreover, since the displace- 
ments are of the same harmonic order, rigid-body motion cannot exist. 

Since the form of the [ K] ^ matrix depends upon the topology 
of the regions , there is some leeway as to the distribution of the zero 
terms within this matrix. This may be accomplished by utilizing various 
numbering techniques for the regions of the structure. Several tech- 
niques may be employed to form the total structure stiffness matrix 
rendering it amenable to facile operation. One technique is to form a 
banded matrix. The numbering for topology of a typical common bulkhead 
tank, so as to produce a banded matrix, as shown in Figure 18. Operations 
with a banded matrix are more efficient, consume less time and less com- 
puter core storage. Another technique presented in Reference 135 does not 
employ banding but is also critically dependent upon judicious topology number- 
ing. For cases where neither of the aforementioned techniques are feasible, 
the stiffness matrix'may be compacted. This may be accomplished by 
eliminating the large amount of zeroes in the matrix from the computer 
storage. This technique in addition to minimizing computer core storage, 
may result in operational simplification. 
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Figure 18 Sample Structural Numbering for Diagonalization of Stiffness Matrix 




Reduced Stiffness: The stiffness matrix established previously must now 
be altered in order to take into account the existence of ring reinforcement, 
possible attachment of the shell to other structures, as well as to inhibit 
rigid body motion and satisfy specific support conditions. 

In the case of ring -stiffened shells, the ring reinforcing matrices 
established in Appendix B, must be stacked in the stiffness matrix in ac- 
cordance with the topology of the ring stiffeners. 

If the shell of revolution under consideration is attached to other 
structures, the stiffness matrix should be modified to take into account the 

effect of these structures. For example, if the shell rests on an elastic 

'■!■■■ ‘ ! j 

support, such as a soil foundation or another structure, the stiffness matrix 

3 : ■ ■ » ' * 

of that support can be written as 


[ F ] « [ K ] [A ] 
m m m 


( 3 - 25 ) 


where sm denotes the joint of the shell at the elastic support. The stiff - 

r ' .. , i ■ i 

ness [K 1 should then be stacked into the total shell stiffness matrix 
m J 

/s « i 

[K],j, at the location corresponding to joint m. The aforementioned tech- 
nique may also be employed in solving problems associated with non- 
axisymmetric structures [ill] being connected to the shell of revolution. 

In the actual shell structure the displacements and rotations of 
the joints of the segments into which the shell must be subdivided, may 
assume specified values, or may be constrained externally. The number 
of displacements which are not specified will be referred to as the degrees 
of freedom of the shell structure. The total number of displacement com- 


ponents specified at the various joints in the structure, will be referred 
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to as the number of boundary conditions. In order to alter (reduce) the total 
the stiffness matrix to take into account the effect of the boundary conditions, 
a boundary condition matrix [BC] must be established. The formulation of 
the matrix may be illustrated by referring to the shell of revolution shown 
in Fig. 19* 

For example, consider the shell of revolution shown in Figure 19a. 
The geometry of this shell suggests the subdivision of the shell into the 4 
segments shown in Fig. 19b. Notice, that the segment between joints 
2-3 may be considered as an equivalent ring stiffener, and its stiffness 
matrix may be computed, and stacked into the total stiffness matrix of 
the shell according to its topology, discussed previously. However, the 
length of the segment 2-3 may be taken as small as desired, whereas the 
lengths of the adjacent segments must increase appropriately to close the 
gap. In the limit, when the length becomes very small, the stiffness matrix 
of segment 2-3 will become a null 'matrix. 

Such a segment will be referred to as a kinematic link. This link 
will affect only the boundary conditions of the adjacent segments. The use 
of a kinematic link where permissible, in lieu of an equivalent ring stiffener 
of chosen finite dimensions, will eliminate the need for computing the stiff- 
ness matrix of this ring. In the example structure, the [BC] matrix in this 
area will reflect only the kinematic relations between joints 2 and 3. 

From Fig. 19b, it is evident that joint 1 is connected to a heavy 

boss. Thus, we may assume that this joint may move solely in the Z and 0 

(tangential) directions. As an alternative, it may be assumed the joint 1 
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is connected to an elastic support. In this case, the stiffness characteristics 
of the support may be appropriately inserted in the total shell stiffness matrix. 
Then, in the formulation of the boundary condition matrix, joint 1 should be ; 
considered totally unrestrained. 

Referring to Fig. 19b, it is apparent that joint 4 is completely 
unrestrained. This joint is merely a point wherein the shell geometry changes. 
Therefore, the [BC] matrix will not impose any constraint on this joint. 

Joint 5 (Fig. 19 b) is provided with an external membrane support. 
Inasmuch as the stress resultants and displacements at the joints of the 
structure are specified in global coordinates, the [BC] matrix will contain 
a trigonometric coordinate transformation for joint 5. 

Thus, the total displacement matrix may be expressed in terms of the 
matrix containing only the actually unspecified (unknown) displacements, 

[£] t = [BC] [£] f (3-26) 

For the example structure of Figure I 9 , Equation (3-26) is given by: 


10b 



yi> 

yi) 

yi> 

V 1 ' 

4 t (2) 

V 2) 

V 2 > 

V 2) 

i T (3) 

A Z (3) 

V 3 > 

V 3) 









The following items should be noted relative to Equation (3-27) • there is a blank 
row in the M matrix for each displacement component or rotation specified as zero 
(fixed) . There are no components for the dependent joint 3 In the matrix [A]p. The 
kinematic relations for this joint are given in the M matrix. The meridional 
component A^ ( 5 ) does not appear since it is fixed, but the perpendicular component 
A^ ( 5 ) contributes to both A z ( 5 ) and A^ ( 5 ). 

By using the definition of work, it may tfe shown that the stress resultants at 
the joints in the directions of the unconstrained displacement components may be 
expressed in terms of the total stress resultants at the joints. This relationship 
is 

[F] f =[BC] T [7'] t (3-28) 

Substituting Equations (3-24) and (3-26) into Equation (3-28), we obtain 


[ 7 ] f = [ BC ] T [ K ] T [BC] U] F + [BC] T [L] t 
R ewriting Equation (3-29) we have 

[F] F = [K] F + [7] f 

where 

[K] f = [BC] T M t [BC] 

W F = [Bcf [L] t 

Inverting Equation ( 3 - 30 ), we get 

t^F = ft F (ff] p - W F ) 

where 


(3-29) 

(3-30) 

(3-31a) 
(3 -31b) 

(3-32) 

(3-33) 


Thus, the total displacement solution for the joints of the structure may 
be obtained by using Equations (3-32 ) and (3-26). 



(3-34) 


Thus for the region ends, combining Equations (3-26) and (3-32) 

, fM = ,[ BC ] (Ml% - 1%) 

an<3. in the interior of each region, for each segment 

[tCr 1 (cCj (m + h)) 

Final Stress Distribution: As noted above, subsequent to obtaining the 
end displacement at any segment, we convert to local coordinates, using 
Equation (3-3) 

[6(i)] = [lDT] T [A(i)3 (3-3') 


The stress resultants at every segment-edge are established in local 
coordinates by combining Equations (3-1) , (3-9) , and (3-11) 


[f (i )] = [IFT f 



A(if 

• • « • 

4j) 


+ U(i)] 


(3-36 ) 


The stress resultants in any elastic restraints may be established from 
Equation (3-25). 

Subsequent to obtaining the stress resultant and displacement 
distribution at the edges of all the segments of the structure, the stress 
resultant and the displacement distribution within each segment must be 
established. This is necessary, inasmuch as in a shell structure having 
a complex shell geometry, the maximum values of the stress resultants 
and/or displacements may not occur at the edges of the segments. Finally, 
it is essential to verify that the established solution satisfies the continuity 
conditions at all the joints of the segments, thus insuring that the errors 

induced during the integration process did not accumulate. 
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The establishment of the stress resultant and displacement dis- 
tribution throughout each segment of the shell, and the checking of the es- 
tablished solutions, may be accomplished simultaneously by an integration 
throughout all the segments of the shell, as described previously in this 
chapter. This final integration, however, does not use the unit vectors 
described earlier as the initial conditions at joint i , but rather the stress 

resultant and displacement vectors obtained from Equations (3-3 r ) and (3-3 6). 

til 

From this integration the stress resultants and displacements of the j 

end are obtained, as well as their distribution throughout the segments. 

The accuracy of the solution obtaiiied by the numerical integration may be 

th 

established by comparing the stress resultants and displacemnts at the j 

' " - til''' 

end of every segment with their counterparts at the i end (same structural 

point), of the corresponding adjacent segments. 

It should be noted that for nonlinear problems the method of 
analysis presented in this chapter must be repeated several times for every 
load increment, as outlined in Chapter 2. After each trial solution the non- 
linear terms in the N e wton - Raphs on procedure are reevaluated using values 
obtained in the previous trial and a check for convergence at the current 
load level is made. Then the load can be incremented again and the pro- 
cedure repeated (see Chapter 2). 

While the current formulation is strictly valid only for shells of 
revolution. Reference 136 has shown how the concepts involved in this 
formulation might be extended to obtain approximate analyses of non- 
circular prismatic shells. 
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CHAPTER 4 


CLASSICAL BUCKLING LOADS FOR SHELLS OF REVOLUTION 
SUBJECTED TO STATIC LOADING 

Various methods were presented in the preceding chapters for solv- 
ing linear and nonlinear static problems for shells of revolution subjected to 
symmetric and unsymmetric loadings. In this chapter, methods for estab- 
lishing the classical buckling load of shells of revolution will be presented. 

The classical buckling load is the load required to bring the idealized "per- 
fect" shell to a bifurcation of its equilibrium (prebuckled) state. That is, 
we shall not be directly concerned either with the postbuckling behavior 0 - 59 ] or 

the effect of initial imperfections on the buckling loads and modes fifio], 

A method was presented in Chapter 2, for establishing the maximum 
value of the load wherein the prebuckling deformation of the shell corre- 
sponding to the applied load becomes unstable. Increments of the load were 
applied to the shell, and using Newton's method, the nonlinear response of 
the shell corresponding to each load increment was established. The maxi- 
mum value of the load for which the prebuckling deformation of the shell 
becomes unstable was established as the point at which the solution ceased 
to converge. In addition to the lengthy computer time involved, this technique 
has other disadvantages. For example, in the case of axisymmetric loading, 
only the n=0 axisymmetric buckling load may be established with this technique 
without the use of "load perturbations". As will be discussed subsequently, in most 

cases of unsymmetric loading, the actual buckling mode may be established with the 
method discussed in Chapter 2. The principal difficulty in this case is the 
extensive computer time required for the analysis. 

It should be noted, however, that the above method may be more useful in 
predicting shell imperfection sensitivity. In the case of axisymmetric imperfect- 
ions for a spherical cap, this was demonstrated in Reference 158, wherein only the 
axisymmetric buckling mode was investigated. However, in combination with the 
"load perturbation" technique, this procedure can logically be extened to study 
unsymmetric modes. In order to study the effects of unsymmetric imperfections, a 
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coupled harmonic geometric and load formulation would be necessary. The beauty of 
the procedure is that no specialized imperfection analysis, other than the 
definition of the stress-free geometry for the most significant imperfection, is 
necessary. The procedure is also independent of imperfection magnitude. Given a 
general nonlinear equilibrium program, imperfection sensitivity may be investigated 
by merely adding another geometry to the program library . 

Derivation of the Stability Equations: The stability equations for a shell of 

revolution can be obtained by the energy methods outlined in Reference 100. 
However, in this investigation, the procedure presented in Reference 137 
will be employed. The typical variables in the general nonlinear Equations 
(1-27, 28, 29) presented in Chapter 1 for homogeneous orthotropic shells, 
will be denoted by Y, and will be decomposed into two components 

Y=Y p +Y B . (.4-1) 

Yp represents the value of the variable at the prebuckled equilibrium state. 

Yg represents the change due to the buckling. The variables Y and Yp must 
satisfy the general nonlinear Equations (1-27, 28, 29). Substituting Equation 
(4-1) into Equations (1-27 , 28, 29) a set of nonlinear equations involving Yp 
and Yg are obtained. The set of prebuckling equations may be obtained from 
the above mentioned set by setting Yg = 0. Subtracting the one set from the 
other, and neglecting terms nonlinear in Yg, the following stability equations 
are obtained (For convenience of presentation the subscript B is omitted in 
all terms except for the "load" term. ): 
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Inasmuchas the variables Y and Yp satisfy the given boundary conditions 
at the ends of the shell, the variables Y_ will satisfy homogeneous boundary 

JJ 

conditions. 


If reinforced or laminated shells are to be analyzed, Equation (4-1) 
must be substituted into Equations (1-31) or (1-32) instead of the corre- 
sponding equations in sets (1-28, 29). Thus, for ring-stringer reinforce- 
ment, the following equations must replace the corresponding equations 
in the set (4-2). 
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For laminated shells the following equations must replace the corresponding 
equations in the set (4-2), 
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Prior to presenting solutions to these equations, it must be decided how 
the prebuckled state will be established. In most earlier buckling analyses, 
this state was established on the basis of the linear membrane theory. This 
procedure yields accurate results for some shell geometries under certain 
boundary and loading conditions, and simplifies the analysis greatly. Re- 
cently, with the introduction of automated numerical analyses, the pre- 
buckled state has been established on the basis of the linear bending theory, 
and even nonlinear bending theory. However, there is not sufficient evidence 
for a general conclusion as to when it is necessary to analyze the pre- 
buckled state by a non-linear analysis. In the last few years, a number 
of shell problems have been solved where nonlinear buckling effects have 
been found to be significant, such as in the case of eccentrically merid- 
ionally stiffened spherical caps [29, 163-165], Most of the general shell 
stability computer programs [28, 117, 138] have options for using nonlinear 
bending analysis for the prebuckling state. 

Stability Under Axisymmetric Loading: As in the nonlinear static 
analysis, considered in Chapters 2 and 3, the stability analysis of shells 
of revolution will be different if the loading is axisymmetric or unsymmetric. 

In the case of axisymmetric loading, the terms in Equations (4-2) having a 
P subscript (prebuckling terms), and the load terms in Equations (4-3) are 
zeroth harmonic amplitudes (invariant with 0). However, the terms which 
do nothave a P subscript (buckling terms), mustbe expressed in terms of the 
Fourier series (2-2). It shouldbe noted, thatin the case of buckling under axisym- 
metric loads, the primed and the double-primed harmonic amplitudes, generally 

will be coupled. However, each pair of harmonic amplitudes (the primed 
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and the double primed) will not be coupled with other pairs and, there- 
fore, each pair may be considered separately. 
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As evident by the presence of the double primed terms in the above equations, 
the primed and double primed terms are coupled, A similar set of equations 
may be obtained involving primarily double primed harmonic amplitudes, with 
only a few single primed cross -coupling terms. From Equations (4-6, 7), it is 
apparent that for certain load conditions (axisymmetric non-torsional loads), 
the pairs of harmonic amplitudes (primed and double primed) may be uncoupled. 
Thus, for axisymmetric non-torsional loads, the buckling modes can be de- 
scribed by the half Fourier series expansions in Equations (2-2). 

In establishing the buckling modes of shells subjected to torsional loads, 
in addition to other loads, the set of Equations (4-6, 7) must be solved simul- 
taneously with the set of equations for the double-primed amplitudes, pre- 
viously discussed. Equations (4-6,7) are valid for orthotropic homogeneous 
shells. If reinforced or laminated shells are to be analyzed, Equations (4-4) 
or (4-5), respectively, must be employed instead of the corresponding equa- 
tions in the sets (4-6, 7). 

It is apparent from Equations (4-6, 7) that for shells subjected to axisym- 
metric loading, buckled shapes may correspond to any harmonic (n = 0, 1, 2* • • ), 
Thus, the buckling loads corresponding to different harmonic buckled shapes 

must be established until the critical (minimum) load is determined. 
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Stability Under Unsymmetric Loading; The appropriate equations for 
the case of unsymmetric loading are analogous to the equations of the non- 
linear unsymmetric static problem presented in Chapter 2. Thus, they 
are more complex than Equations (4-6, 7), for the case of symmetric load- 
ing. For unsymmetric loading, when the loads and displacements are expanded 
in Fourier series (see Equations (2-1,2)) the resulting equations involve 
product terms (see Equations (2-6a,b)). Thus, by referring to Chapter 
2, the general equations for the stability of homogeneous orthotropic shells 
of revolution under unsymmetric load may be written as 
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where, as in Equations (2-8), the above equations have been obtained by utili- 
zing only the half of the series (2-1) with primed amplitudes. A similar set of 

equations may be obtained by using the half of the series (2-1) with the double 
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primed amplitudes. In Equations (4-8), the nonlinear terms involve cou- 
pling of the primed and double primed amplitudes, as well as the harmonics. 
The amplitudes of the harmonics in the nonlinear terms with a subscript P 
represent the effect of the prebuckled state. The nonlinear terms in Equations 
(4-8) may be obtained from Equations (2-9) with the appropriate addition 
of thq P subscript. The buckling "load" terms -(fg^^, f^g^, f^ n ^) are 
obtained from Equations (2-9) by omitting the first single term (Fq^, 

F^ n ^). If reinforced or laminated shells are to be analyzed, Equa- 
tions (4-4) or (4-5) should be used in place of their counterparts in Equa- 
tions (4-8). 

Stability Considerations ; In Chapter 2, itwas noted thatthe resulting 
equa tions couldbe significantly simplified if a line of symmetry existed in the 
loading pattern. This was accomplished inasmuchas,, in this case, the resulting 
deformation and stress pattern is symmetric with respect to this line of symmetry. 
In the case of buckling, if a line of symmetry exists in the loading pattern, 
and consequently in the prebuckling deformation and stress state, the load- 
ing pattern and the prebuckling deformation and stress state may be repre- 
sented by only half of the Fourier series expansions. However, this does 
not denote that the buckling deformation shape may be represented by the 
same half-series expansions. Thus, the full series expansions must be 
used resulting in two different sets of equations corresponding to buckling 
modes "in-phase" and "out- of- phase" with the applied load. This was also 
the case for axisymmetric loading (no torsional loads). However, since 
in that case, each harmonic buckling shape could be investigated separately, 
the two sets of equations represent buckling shapes differing only by a rigid 
body motion. Thus, only one set of equations was sufficient for the analysis. 

If the shell was subjected to torsional loads, in addition to other axisymmetric 


121 



loads, the resulting sets of equations will be cross-coupled and must 
be solved simultaneously for each harmonic buckling mode. In the case of 
unsymmetric loading with a line of symmetry, the prebuckling harmonic 
amplitudes (e. g. primed) are coupled with the buckling harmonic amplitudes 
(primed or double -primed). Thus, the two sets of equations yield different 
buckling loads corresponding to in-phase or out-of-phase buckling shapes. 

In Figure 20, the axisymmetric load versus the resulting deformation 
is plotted, whereas the unsymmetric load versus the resulting deformation 
is plotted in Figure 2l. Under axisymmetric loading, the iteration technique 

liii 

described in Chapter 2 would proceed along line OA BCD which corre- 

•. ' t 

sponds to the axisymmetric nonlinear static analysis. Point D is established 
as the buckling load referred to by Thompson [139] as the "snapping load". 
However, the actual (lowest) buckling load may correspond to a non-axisym- 
metric buckling configuration. ( As shown in Fig. 20, point A , obtained on 
the basis of the stability analysis to be presented in this chapter, for the 
n = n^ buckling configuration may correspond to a buckling load lower 

i 

than that corresponding to point D . ) Thus, for any assumed buckling con- 
figuration corresponding to an assumed value n, the solution to be presented 

i i t 

in this chapter will yield a buckling load (points A or B or C ), whereas, 

t 

the method presented in Chapter 2 will yield only point D corresponding to 
n = 0. That is, the solution to be presented in this chapter permits a buck- 
ling configuration described by harmonics different from n = 0 (describing 
the applied loads), whereas, the solution presented in Chapter 2 yields an 

n=0 buckling configuration only, without resorting to "load perturbations". 

From experimental evidence, however the buckling configuration of shells of 

Notice, that in the case of axisymmetric loading, as discussed previously, 
the harmonics uncouple; consequently, the buckling configuration will corre- 
spond only to one harmonic. 
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Figure 20 shell St ah ility-Axi symmetric Load 
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revolution subjected to symmetric loads is not always described by the same 
harmonic as the applied loads. The two methods for buckling analysis to be 
discussed in this chapter nay be employed to establish buckling loads for any 
buckling mode. That is, referring to Fig. 20, to establish points A', B', C', D' 
if the prebuckling state is described by the nonlinear bending theory, or points 
A, B, C, D if the prebuckling state is described by the linear bending theory. 

I 

Under unsymmetric loading, (see Fig. 21) point A may be established 
by the iteration technique described in Chapter 2. Note, that even if the 
applied load was described only by a few Fourier harmonic amplitudes, the 
iteration procedure would eventually involve all the harmonic amplitudes 
of the stress resultants and displacements. In the methods to be described 
in this chapter, if a nonlinear prebuckling state is to be considered, all the 
harmonics couple, and point A is established. If the prebuckled state is 
analyzed by linear bending theory, several buckling configurations described 
by different families of harmonics constitute mathematically acceptable solu- 
tions. The loads corresponding to these configurations are denoted in Fig. 

23. by points A, B, C, D, E, F. These families of harmonics will be described 

I 

subsequently. The buckled configuration corresponding to point A is 
described by all the harmonic terms, whereas, the buckling configurations 
corresponding to points A, B, D are described by groups of harmonics each 
of which is contained in the terms describing the configuration corresponding 

i 

to point A . Thus, the loads corresponding to points A, B, D will be approxi- 
mations of different order of accuracy to the critical buckling load defined 

i 

by point A , rather than buckling loads corresponding to different buckling 
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configurations as in the case of a shell under axi symmetric loading. 

i 

The load corresponding to point A may be obtained by incremental 
Newtonian iteration, using a larger number of harmonic terms for each 
successive load increment. Thus, the analysis is lengthy. Moreover, in 
the analysis, the Fourier series describing the buckling configuration must 
be truncated, consequently the load corresponding to point A may be es- 
tablished approximately. Points A, B, D, however, may be obtained by one 
eigenvalue analysis (for each point). In this analysis, the Fourier series harmonic 
families describing the buckling configurations must be truncated. If a 
dominant harmonic group exists in the description of the actual buckling 
shape, and the same number of Fourier terms is retained in the linear 
and non-linear analyses, the buckling load obtained on the basis of the 
linear analysis (point A) may actually be closer to the true buckling load than 
the load estimated on the basis of the non-linear analysis. Furthermore, 
the computer time required for obtaining the buckling load on the basis of 
the linear analysis is much smaller than that required for obtaining the 
buckling load on the basis of the non-linear analysis. It should be noted, 
that the smallest in-phase buckling load estimate (point A) will be larger 
than the actual in-phase buckling load. However, the smallest out-of-phase 
buckling load estimate (point C) may be smaller than the smallest in-phase 
buckling load estimate (point A) since they are estimates to different possi- 
ble critical loads (in-phase and out-of-phase). 

Solution by the Determinant Evaluation Method: The determinant 
evaluation method represents the classical approach to the stability problem. 
For any assumed value of the applied load, a static analysis is performed to 
establish the prebuckling stress and deformation components at every point 
chosen for Runge Kutta integration in every shell segment. This is ac- 
complished by employing the me shods described in Chapter 3 for either linear 
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or nonlinear theory. The results of this analysis are introduced as the 
terms with subscript P in Equations (4-6) and (4-7) if the load is axisym- 
metric, or Equations (4-8) if the loading is non-axisymmetric. These 
equations are then integrated using the the Runge-Kutta procedure in a 
fashion analogous to that of the static analysis of Chapter 3 and the pre- 
stressed stiffness matrices of the segments are formed. Note, that a load 
matrix associated with the stability analysis does not exist, inasmuch as 
the loads have been eliminated from the pertinent equations by taking into 
account that prior to buckling, the shell is in a state of equilibrium under 
the influence of the buckling loads. That is, the buckling state is a possible 
second state of equilibrium under the buckling loads. The prestressed , 
stiffness matrices of each segment are 8 x 8 matrices for shells subjected 
to axisymmetric loading, or 8M x 8M matrices for shells subjected to un- 
symmetric loading, where M is the number of harmonic terms retained in 
the analysis. 

The prestressed stiffness matrices of the segments are stacked to 
obtain the total matrix of the structure by a procedure identical to that 
discussed in detail in Chapter 3. This matrix is then reduced by employ- 
ing the boundary conditions. Thus, we obtain an equation analogous to 
Equation (3-30) in the form 

[Kj [A] = 0 (4-9) 

f F F 

or 

detpT' ] =0 . (4-10) 

P F 

If the assumed load was the correct buckling load corresponding to the 
assumed harmonic buckling configuration. Equation (4-9) would be identically 
satisfied. Otherwise the determinant would not vanish. 
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Inasmuchas the pre stressed stiffness matrix may be large, the evalua- 
tion of its determinant on a computer may involve overflow or underflow. 

To avoid this problem, the matrix [K^] is first converted into an upper 
triangular matrix by a technique such as Gaussian elimination [118]. 

The value of the determinant of such a matrix is the product of its diagonal 
terms. Since the determinant vanishes, these terms can be normalized 
by dividing each term by its absolute value. Thus, the value of the determin- 
ant is always ±1. A sign change of the determinant between two consecutive 
loads signifies that the value of the determinant vanishes between these two 
loads. This technique avoids the establishment of spurious sign changes 
[25], 


The assumed load is incremented until the determinant changes sign. 

The load increments may be either constant or varying. The latter may 
be established by extrapolation from the previous load increments. 

When the prebuckled state is analyzed by the linear bending theory, 
only one static solution for one load is required. The prebuckling stress 
resultants and displacements corresponding to other values of the applied 
load may be established by superposition. When the prebuckled state is 
analyzed by the nonlinear theory, static solutions are necessary for each 
assumed value of the load. These solutions are established by the Newtonian 
method described in Chapter 2. Aside from the additional accuracy of a 
nonlinear prebuckling solution [29] , other flexibility of analysis trey be 
gained by including such an option. For example, while normal solution by a 
nonlinear prebuckling analysis does not include consideration of the local 
panel problem considered by Dickson and Brolliar [141] , the similarity 
between the iterative natures of both solutions make a combination possible. 
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Thus, the consideration of the local panel stability problem prior to overall 
shell buckling, can be incorporated into the iterations needed for a nonlinear 

. . . ... . t .... I . \ , ,. ■- , -;L 

prebuckling analysis, ! ; 

With the aforementioned procedure, all the possible buckling configura- 
tion may be checked. In the case of axi symmetric loading, each harmonic 
can be checked independently of the others. Since all harmonics (n = 0 - ~) 
cannot be investigated, the problem remains to insure that the lowest 
buckling load is associated with one of the harmonics checked. It should 
be indicated that an automated checking procedure may result in erroneous 
conclusions. Figure 22 presents the lowest buckling loads, obtained for 
different buckling configurations corresponding to the indicated harmonics 
for the classical problem of a circular cylinder subjected to end compression 
[142], If an automatic procedure is programmed to establish a relative mini- 
mum within a given range of harmonics, the buckling loads corresponding 
to any of the harmonics n = 2, 7, 9, 11 could be obtained as a solution, where- 
as, the actual buckling load corresponds to n = 2. In this example, the values 
of the buckling loads obtained for n = 2, 7, 9, 11 do not differ appreciably, 
however, for each harmonic the buckled configurations differ considerably. 

A basic difference in the stability analysis ofshells subjected to axisym- 
metric and unsymmetric loading pertains to the type of buckling shapes that 
must be considered. If the nonlinear Equations (2-9) were converted to a 
form suitable for stability analysis, as previously discussed in this Chapter, 
and if only one load harmonic (n 4 0) is considered, several buckling con- 
figurations described by different families of harmonics are found as mathe- 
matically acceptable solutions. The prestressed stiffness matrices corre- 
sponding to these different buckling configurations assume the form of the 

129 






0123 4.5678 



Case I: Axi symmetric loading (& = 0) [all 

harmonics are uncoupled in the 
pre stressed stiffness matrix] 
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Case II: Antisymmetric loading (& = 0,1 ) 

[all harmonics are coupled in the 
pre stressed stiffness matrix] 
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Case III: Unsymmetric loading (i = 0,2) [Prestressed stiffness matrix reduces to two 

Uncoupled families of harmonies] 
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Case IV: Unsymmetric loading (£ = 0,3) [Prestressed stiffness matrix reduces to two 

uncoupled families of harmonics] 
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Case V: Unsymmetric loading (Z = 0,4) [Prestressed stiffness matrix reduces to 

three uncoupled families of harmonics] 


Figure 23a Forms of the Prestressed Stiffness Matrices Corresponding to 
Single Unsymmetric Harmonic Loads 
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Figure 23b Forms of the Prestressed Stiffness Matrix Corresponding to Three Harmonic 
(0,2,4) Loading. 
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trident matrices, shown in Fig. 23a, b. Figure 23 shows how these pre- 
stressed stiffness matrices corresponding to different harmonic loadings 
uncouple into families of harmonics. Each family of harmonics represents 
a different estimate of the buckling load (see Fig. 21, points A, B, G, D* • • ), 
each associated with a different buckling configuration. If one of these con- 
figurations is a close approximation to the actual buckling configuration, 
then the buckling load corresponding to this configuration will be lower than 
all the other estimates (points B, C, D* • • in Fig. 21) obtained on the basis 
of the linear prebuckling analysis. Furthermore, this estimate may be lower 
than the buckling load obtained on the basis of a nonlinear prebuckling analysis 

(point A*, Fig. 21) retaining the same number of terms in the Fourier series. 

In Reference 118 a stability analysis of a shell of revolution sub- 

1 

jected to some types of unsymmetric loading is presented. However, in 
this reference several nonlinear terms usually included [20, 24, 28, 32, 117] 

in shell buckling analyses have been omitted. Also, only the family of 
harmonics which are multiples of the load is assumed to represent the 

configuration corresponding to the lowest buckling load. This assumption is 

based on the conclusion (invoking St. Venant's principle [118, p. 77]) that the 

effect of load harmonics with -t >2 dissappear at a small distance (within a 

diameter of the latitude circle) away from the loaded edge. This conclusion, 

however, is not always valid. Transverse loads, such as have a 

short decay length; however, the decay lengths of the in-plane loads such 

( 2 ) ,( 2 ) 

as ' or T^q may be many times larger than the diameter of the latitude 

circle, even for shells with ring reinforcement [111, 143]. In certain cases, 

the assumptions of Reference 118 may lead to contradictory 

conclusions. For instance, consider a shell which buckles in the n = 10 

harmonic configuration when subjected to axisymmetric load. If an J = 3 
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harmonic load is added to this shell, the assumptions of Reference 118 
lead to two possible conclusions. If it is assumed that the additional load 
(t= 3 harmonic) represents only an edge effect, then the buckling configura- 
tion will be similar to that (n = 10) established when the shell was subjected 
to axisymmetric loading. Thus, the best estimate for the buckling load tn^y 
correspond to a buckling configuration associated with the family of harmonics 
which contains n= 10 (see Figure 23a, b, Case IV). On the other hand, if the 
buckled configuration corresponding to the lowest load is obtained from the 
family of harmonics which are multiples of the loads, then the buckled con- 
figuration corresponding to the lowest buckling load in the aforementioned 
example should be obtained by n = 0, 3, 6, 9, 12, 15. . . . 

The added complexity of the unsymmetric stability analysis, over the 
axisymmetric, can be seen from the above discussion. Even though only 
one load harmonic (l) is applied, and the prebuckling analysis is on the 
basis of linear bending theory, families of harmonics must be employed to 
establish the estimates for the buckled configurations. Thus, larger ma- 
trices are involved, and must be treated as discussed in Chapter 3 for 
the nonlinear unsymmetric static analysis. If the nonlinear theory is em- 
ployed in the prebuckling analysis, all harmonics must be used to describe 
the buckled configuration. 

The number of terms to be retained in the truncated Fourier series ex- 
pansions of the shell functions requires further investigation in both cases 
of stability analysis of shells subjected to unsymmetric loading, and cases of 
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unsymmetric nonlinear stress analysis using the Newtonian iteration tech- 
“ niques [20, 22]. It may be possible to obtain accurate estimates of the buck- 
ling configuration by retaining only the n = 0 harmonic and a few more judi- 
ciously chosen harmonics from a family of harmonics, not necessarily in 
"a consecutive order. If this were the case, itwouldbe evenmore preferable to employ 
the eigenvalue approach rather than the Newtonian iteration procedure. In 
this case, a few terms from the dominant harmonic family would yield more 
satisfactory results than the same number of harmonic terms in the New- 
tonian iteration solution, since more significant terms would be contained in 
the eigenvalue solution. For instance, some shells such as relatively shal- 
low ellipsoidal heads subjected to internal pressure, buckle in a relatively 
high harmonic pattern (n~ 50). If the ellipsoid is subjected to an unsymmetric 
load, in addition to internal pressure, the question arises as to whether the 
buckled configuration can be approximated by some lower harmonics and 
some in the proximity of n = 50, omitting the intermediate members of the 

harmonic families, or all the harmonics up to n= 50 must be retained for a 
satisfactory approximation of the buckled configuration. 

It should also be noted, that the high local wrinkling associated with high 

harmonics, may require the inclusion of shear deformation in the theory 
employed [113]. 
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In conclusion, it maybe stated that the determinant evaluation method, 
generally, has two basic disadvantages. Firstly, very* extensive computer 
time is required for establishing the buckling load of shells of complex 
geometry. This is primarily due to the fact that if the magnitude of the 
buckling load cannot be estimated apriori^ many load increments may have 
to be considered before the vajue of the load causing the determinant of 
the prestressed stiffness matrix to vanish is established. The second prob- 
lem is even more complex, since it is possible that the low buckling loads 
are close together [145] even for buckling configurations described by a 
single harmonic. Thus, a load increment may skip two close roots without 
causing the sign of the determinant to change. To circumvent the foregoing 
disadvantages, other stability analysis techniques have been developed. 

Eigenvalue Methods: When a problem involving the stability analysis 

of any structure is analyzed by finite element methods, it may be reduced to 
a linear eigenvalue problem of the form 

([A] + X[B]){a] = 0 (4-11) 

where [A] is the stiffness matrix of the elements and [B] is the incremental 
stiffness matrix. Each of these matrices is formed separately by assuming 
a displacement function for the element. In finite difference or numerical 
integration techniques, however, a linear eigenvalue formulation is not 
readily deduced. The first eigenvalue -type analysis employing numerical 
integration techniques, was formulated by Cohen for natural vibrations of 
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shells [113], and subsequently for stability analysis of shells subjected to 
axi symmetric loads [117]. A similar method was later utilized in Reference 
31 for finite differences. The method is iterative, based on the Stodola 
technique [114], and is essentially the inverse power method [11 5]. 
Basically, the homogeneous equations resulting from stability (vibrations) 
analysis are converted into a series of nonhomogeneous equations by assum- 
ing a buckled (vibration) shape and, thus, creating nonhomogeneous terms. 

The solution of this problem provides a more satisfactory estimate for the 
buckled configuration which, in turn, is employed to establish a new set of 
non- homogeneous terms. The procedure is repeated until the lowest eigen- 
value (buckling load or frequency) is obtained for the harmonic configuration 
under consideration. 

This method requires less computer time than the determinant evaluation 
method, and moreover, the possibility of skipping roots is eliminated. How- 
ever, it has some disadvantages. In vibration problems, wherein it is 
necessary to compute higher frequencies, in order to establish the higher 
roots, all the lower roots must be swept out [114]. It has also been found, that 
the time required to establish two consecutive roots is a function of the ratio 
of the value of the lower root to the higher root. That is, the time is larger 
when this ratio approaches unity. Thus, in order to decrease the time re- 
quired for establishing the second root, the origin should be shifted to the first 
root [116]. In stability problems, where only the lowest buckling load for a 
particular circumferential configuration is of interest, the aforementioned 
drawbacks of the inverse power method of solution do not exist. However, 
for a shell of complex geometry, convergence may be slow depending upon 

the initial choice of the eigenvectors (u, v, w) [114, 140]. 
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In a recent investigation [138], the finite difference technique is applied 
to the shell energy integral, rather than the differential equations of equilibri 
um and separate stiffness and incremental stiffness matrices are formed. 
However, the inverse power method is still employed to solve the resulting 
linear eigenvalue problem. 

To overcome the foregoing difficulties, a different formulation of C j 
stability problem has been presented in Reference 157. The basic pre- 
stressed stiffness matrix, [K ] of the structure, is an unknown transcen- 

P F 

dental function of the prestress state variables. In vibration analysis, 
Przemieniecki [146] shows that the dynamic stiffness matrix is actually an 
infinite power series on the frequency. This finding may be extended to 
stability analysis. Thus, we obtain 


[iT] = [k] + xtf£] 

P F F 1 ] 


x 2 [k 


II J 


+ \ 3 [k 


III 


] + 


(4-12) 


where X is the buckling load. In Reference 146, it is shown that the ratio 

of consecutive matrices [K.] / [K. , ,] is of the order of Young's Modulus 

1 F 1+ F 

of the structure. In formulating the stability problem by finite element 

methods (see Equation (4-11)), the matrix [A] is an approximation to [k] f , 

whereas the matrix [B] is an approximation to [K-] , to the order of accuracy 

1 F 

of the assumed deflection functions for the element. In the numerical inte- 
gration technique £k]_, and [K ] can be formed exactly using the exact 

* P f 

differential equations. However, it is impractical to form the other matrices. 
Inasmuchas the relative magnitudes of these matrices are known, the follow- 
ing solution technqiue has been proposed [157]. Using Equation (4-12), the 
stability Equation (4-9) can be cast into Hie following form: 

{[K] F + \[kJ] + X 2 DK^] + \ 3 [kJJj] + • • • 3 [A] F = 0 

F F F 
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(4- 13a) 




or 

t(K] F + 3 ^-([K^T^')] - ® f )}[4] p = o (4-14) 

where X^ ^ is the buckling load estimated in the (i - l)st trial. 

The iteration Equation (4-14) is utilized as follows. As in the deter- 
minant evaluation method, the prebuckling analysis for the establishment 
of the prebuckling terms in Equations (4-6) or (4-8) is performed for a 
chosen value of the load on the basis of either linear or nonlinear bending 
theory. As in the determinant evaluation technique, the reduced prestressed 

is formed, where in this no- 
F 

tation, X^ ^ signifies the chosen value of the load. The structure stiffness 
matrix, [K]p (without the prestress terms), is also formed, for the buck- 
ling configuration under consideration. The subtraction of these two ma- 
trices, as in Equation (4-14), isolates that part of the prestressed stiffness 
matrix which is dependent upon the buckling load. A linear eigenvalue 
(Equation (4-14)) problem is then formed, and solved for the new value of 
the load, X^. The iteration sequence converges when X^ = X. ^ or X^/Xj_j 
approaches unity to a desired degree of accuracy. Although this must be 
accomplished for each root desired, satisfactory approximations to higher 
roots are also available when X^/X^ ^ reaches unity for one root. This is 
due to the fact that the numerical integration technique does not lead to large 
matrix equations and, thus, there is computer storage available for an eigen- 
value solution algorithm which provides all the roots of Equation (4-14). In 
this case, the convenient in-core algorithm used is the Householder [115] 
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stiffness matrix of the structure [K (X^ j)] 



method. Cohen [117] has shown that the numerical integration formulation 
of the stability problem produces real roots. This indicates that the ma- 
trices in the linear eigenvalue problem will both be symmetric, and one 
will be positive definite. Thus, the applicability of the Householder tech- 
nique is assured [115]. 

The specialization required for the consideration of unsymmetric 
loading as opposed to axisymmetric loading, and nonlinear prebuckling 
analysis as opposed to linear prebuckling analysis, is applicable to this 
method as well as to the determinant evaluation procedure discussed pre- 

.'■’S 

viously. In this case, it should be noted that the stiffness matrix [K]^ 

will not involve coupling of the harmonics, regardless if the load is symmetric 

or not. The [K ] matrix, however, will be coupled as described pre- 
P F 

viously. 

An iteration equation analogous to Equation (4-14) has been presented 
by Bushnell [138] for establishing the critical load at buckling of shells of 
revolution, using a nonlinear prebuckling state and finite differences. This 
Equation of Reference 138 may be obtained by rewriting equation (4- 13b) 
as, 

{[KjT^)] +\.([K^r^)] - [K] f )3 [A] f = 0 . (4-15) 

F F 

In this case, for convergence X, -* 0. It should be noted, that in the formu- 
lation herein presented, iteration is necessary regardless of whether the 
prebuckling analysis is linear or nonlinear. This is due to the fact that the 
prestressed stiffness matrix is defined by Equation (4-12), whereas, in 
Reference 138 for the case of a linear prebuckling analysis, Equations (4-15) 
reduce to Equations (4-11), since the prestressed stiffness matrix is 
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approximated as being linear in the eigenvalue. Either of Equations (4-14) 
or (4-15) can be employed to obtain a solution to the stability problem, after 
few iteration cycles. Examples using Equations (4-14) will be presented in 
Chapter 6. 

Many stability problems have a double -loading system. For example, 

if a tank with an insulating wall is manufactured at room temperature and 

then partly filled with a cryogen, the tank is subjected to a state of stress. 

If this tank is then accelerated, it is subjected to mechanical loads which 

may cause buckling. The effect of the thermal prestress can easily be 

considered in the analysis, by including it in the [K]p and [K (X. _ ^ )] ma- 

P F 

trices of Equations (4-14). This will necessitate the solution of two static 
prestress problems, one with the thermal loads and the other with com- 
bined thermal and mechanical loading. Therefore, even for a linear pre- 

\ 

buckling analysis, for shells subjected to unsymmetric load, the harmonics 
will couple in the matrix [K]p as adjusted for the thermal effect. 
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CHAPTER 5 


NATURAL VIBRATIONS 

In this chapter, the free vibrations of shells of revolutions from, a 
stress -free or a pre stressed state are analyzed. The results of this 
analysis may be applied to establish the dynamic response of shells of 
revolution subjected to a harmonic exiting force or to any transient load- 
ing if the modal approach is to be employed. 

Dynamic Equilibrium: The nonlinear equations of motion of shells 

of revolution may be obtained from the equilibrium Equations (1-20) by 
including the effects of meridional, circumferential, normal, and rotatory 
inertia. Thus, we obtain 

r l N 0,9 + ir (N cp0 r o 2) 'cp- Q 0 r l sin ^= - r l r o (f 0 V + f 0 ) + r l r o (a O U ' a lV 
( Vo K cp + r l N cp0,0 - N 0 r l cosCp - r o Q cp = ■ r l r o (f <^ + V + r l r o (a O^ + a l S 0 ) 
( Q cp r o ) 'cp + r l Q 0,0 + r o V N 0 r l Sincp= " r l r o (f C* + V + r l r oV" (5_1) 


- r l M cp0,0 - ^cp^'cp + M 0 r l cosCp+ r l r o Q cp = " r l r o m 0 + r l r o (a r + a 2 a 0 ) 

”^ M cp0 r o ), cp - r l M 0,0- M cp0 r l cosC P + ^r o Q 6 = "Vo" 1 ** r i r o (a l U “ a 2V 

where the (* ) signifies differentiation with respect to time; the load 

L(i = 0,cp, Q) and nonlinear terms £. (i = 0,cp, Q) are defined in Equations (1-21) 

and (1-22) respectively, and 


i. = J pC (j) dC 


(5-2) 


p(C) being the mass density. It is recognized, that generally, if the effect 

of rotatory inertia is significant, the effect of shear deformation is not 
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negligible. However, in this analysis, as in References 113 and 118, the 
effect of shear deformation will be neglected. In reinforced shells, the 
combined rotatory inertia of the skin and reinforcements may influence the 
results even when the effect of shear deformation is negligible. For mono- 
coque or sandwich shells, a^ = 0 since in these shells, C is measured from 
the centroid of their cross-sections. For reinforced or laminated shells, 
the a. will all be different from. zero. 

J 

Utilizing Equations (5-1) and following a procedure analogous to that 
presented in Chapter 1, a set of equations analogous to Equations (1-27) but 
including the inertia effects may be obtained. If a reinforced or a laminated 
shell is to be analyzed, Equations (1-31) or (1-32) may be employed without 
any modification. As in Chapter 4, in the case of stability of shells, the 
typical shell function denoted by Y will be considered as the sum of its 
values Yp, in the prestressed equilibrium state and its change due to the 
vibrations, Yye lU,t . Thus, we have 

Y = Y p + Y y e m (5-3) 

where i = and 0) is the circular frequency of vibration. 

We shall postulate that the rotations due to the vibrations are small 
as compared to unity, and consequently, we shall disregard the terms in- 
volving products of these rotations with the stress resultant amplitudes due 
to the vibrations. Thus, following a procedure analogous to that described 
in Chapter 4, we obtain: 
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These equations can be integrated numerically to establish the frequencies 
and mode shapes of free vibration of shells of revolution subjected to axisym- 
metric or unsymmetric prestress. The numerical integration method has 
been employed only in establishing the frequencies of non-prestres sed shells 
of revolution [112, 113]. A finite difference solution for axisymmetric states 
of prestress has been presented in Reference 24. 

In the analysis of a reinforced or laminated shell, Equations (4-4) or 
(4- 5) may be employed without any modification. 

In the above equations, the term e lWt has been factored out. Thus, 
the solution of these equations will yield the amplitudes of stress resultants 
and displacements. It should be noted, that in Equations (5-4, 5) the terms 
involving products of rotations due to the prestress and stress resultant 
amplitudes due to vibration are retained. It is realized, however, that the 
rotations due to the prestress are small as compared to unity, and conse- 
quently their products with the stress resultant amplitudes due to vibration 
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may be negligible as compared to the stress resultants. 

Vibration Under Prestress : For axisymmetric prestress, the pre- 

stress terms are of the zeroth harmonic, while the vibrations may involve 
any single harmonic of the Fourier components. For non -axisymmetric pre- 
stress, the prestress harmonics will couple with the vibration state harmonics 
in a fashion analogous to that discussed in Chapter 4 for buckling of shells 
subjected to non- axisymmetric loads. Inasmuchas the formulation of the 
problem of stability and the problem of vibration under prestress is similar, 
the analysis of the special cases of prestress presented in Chapter 4 is 
valid for the analysis of vibrations under prestress. Hence, Equations (4-6) 
may be employed for problems involving vibration under axisymmetric pre- 
stress, whereas, Equations (4-8) maybe employed for problems associated 
with vibrations under non- axisymmetric prestress. These equations must 
be modified, firstly, by adding to the terms f.g (i= 9 , cp, Q) the effects of in- 
ertia. Thus 


f (n) 

± 0V 


f (n) 

cpV 


f (n) 

*CV 


o 

' £ t pB <n)+a,2<a O V<n)+a i n 0 (n)) 

- f $B (n) + ll) 2 (a 0 wt n > - ^-[najU (n> - . 


(5-6) 


Secondly, by adding the following term to Equation (4-6d) for vibrations under 
axisymmetric prestress and to Equation (4-8d) for vibrations under non- 
axisymmetric prestress 

+ou 2 (a 1 V (n) + a 2 n 0 (n) ) . (5-7) 
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The variables with a subscript B are defined in Equatipns (4-7) for vibrations 
under axisymmetric prestress aiid are obtained from Equations (2-9) for vi- 
brations under unsymmetric pre stress. 

If a reinforced or a laminated shell is to be analyzed, Equations (4-4) 
or (4-5) may still be utilized without any modification. 

For free vibrations of non-prestressed shells all the terms with a 
subscript P, as well as the prestress loads F. (i = 9, Cp, Q) vanish. 

Numerical Solutions : The same two methods used for solving stability 
problems will be used in solving vibration problems. Thus in the sequel only 
a brief discussion of the application of these two methods will be presented. 

Determinant Evaluation Method: In problems of vibration about a 

prestressed state, the necessary prestress terms must first be determined. 
This is done by means of the static analysis outlined in Chapter 3, using 
either linear or nonlinear theory. A value of the frequency is then assumed, 
and the vibration Equations (5-4,5) are utilized to form a stiffness matrix 
as in Chapter 4. In the case of vibrations this "dynamic stiffness matrix" 

[146] is a function of frequency. The natural frequencies are established as 
the frequencies which render the determinant of the dynamic stiffness matrix 
for the structure equal to zero. The techniques of finding these frequencies 
are the same as for finding the critical loads at buckling. 

Linear Eigenvalue Methods: The two techniques discussed in Chapter 

4 for stability are both applicable to the vibration problem. The equations, 
corresponding to Equations (4-14), suitable for analyzing vibrations under pre- 
stress can be formulated in the following manner. The prestress stiffness 

matrix, [K ] , of the structure, is formed from the static solution for the 

P F 

shell subjected to the given pre stress. The dynamic stiffness matrix, [K_.] , 

F 

ibj 



of the structure, is formulated by assuming a value for the frequency and 
using the dynamic Equations (5-4, 5). Thus, following the stability formu- 
lation as discussed in Chapter 4, we may write 


00. 

i_ 

pJ F uu? 


CDKJ + 




i-1 


F 


- [Kj )UA] ir = 0 


(5-8) 


These equations may be solved by the eigenvalue solution technique and the 
iteration procedure discussed in Chapter 4. It has been shown [113] that 
the numerical integration method for solving the vibration problems of 
shells of revolution yields real, positive frequencies. This indicates that 

I 

in Equation (5-8) both matrices will be symmetric and positive definite. 
Thus, the applicability of the Householder technique is assured [115]. 


It should be noted that Equation (5-8) could be reformulated as follows, 


i>i p + < i )] F - [ s ] F ] 5 [a] f = 0 • 


(5-9) 


In this case, convergence is indicated as tt). -* 0 while in Equation (5-8) con- 
vergence is indicated ae ot)? / ou? ^-*1. 

The advantages of using the Householder technique and of using the 
formulations (5-8, 9) over the determinant method were noted in Chapter 4. 

Of specific importance in dynamics of shells is the ability to quickly estimate 
many frequencies while having completed an iteration solution for only one. 


For problems of free vibrations of unpre stressed shells of revolution 
the static stiffness matrix is used in place of the prestressed stiffness matrix 
of the structure. All other operations remain the same. 



Critical Speeds of Rotating Shells : The matrix equations for shell 


dynamics, referred to a coordinate system which rotates with the shell about 
its axis of revolution with constant velocity, ft, may be cast in the form 

+ 2Q [d] f {A'k +(£], - & ={o} (5-10) 


where {a*^ and jKp J ^ have been defined previously as the incremental displace- 
ment vector, and the fixed preload stiffness matrix respectively, and where 

A 

C^pd 3 j, is the variable load prestress stiffness matrix including the rotating 

A 

inertia terms of the incremental state. The additional matrices, [Ml 

A 

f D } j, contain the nonrotating inertia and the Coriolis dynamic effects, 
respectively. The dynamic stability of a rotating shell may be investigated by 
substituting 



F 



(5-11) 


into Equation (5-10) to obtain 

-U 2 [M] r {A} p + 21(0 o [D] p {a} f +( D^J r = {0} (5-12) 


In the presently programmed version of the STARS vibrations program ( STARS -2V ) , 
Equation (5-12) is not Solved in complete form but rather for special cases. 

For the vibration case discussed previously in this chapter, ft = 0, and Equation 
(5-12) will reduce to Equation (5-8). The code words FREV and VPRE will then 
define the shell as being either stress-free or axisyrametrically prestressed. 

For the case of critical speed analysis, to= 0, and Equation (5-12) 
reduces to 

< &1f - ^ > Mf - {«} 


(5-13) 



In this case the matrix [K^] p contains the effects of the centrifugal 
accelerations on both the prestressed and incremental deformation states of 
the shell. The code words CRSP and PCRS again are used to define the shell 
as being initially stress-free or axisymmetrically prestressed (static load). 

At the request of NASA fl66l a third option is also available wherein 
00 = fQ. In this case the Coriolis effects are neglected and the Equation 
(5-12 ) reduces to 

-M [m] f {4} f + ( [Bp], - Cp [k^)],) {i} r = { 0 } (5-14) 

The neglect of the Coriolis terms is not a serious violation and is consistent 
with the fact that torsional prestress in' not allowed in the present programs. 

In addition, the effect of the Coriolis terms has been found to decrease for 
higher rotational speeds [ 16 T] such as those of interest for the present 
analysis. The above option may be utilized to approximate cases in the analysis 
of rotating shells wherein some additional dynamic load of unknown description 
is causing a "mass perturbation" upon the shell. This perturbation is described 
proportionally through the multiplier upon the critical speed, f, and is applied 
in whatever harmonic is being investigated for critical speeds. Again the code 
words CRSR and PC SR differentiate only between the non-existence or existence 
of static preload. 
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CHAPTER 6 


NUMERICAL EXAMPLES 


In this chapter, solutions using the programs herein documented will he 
compared to solutions utilizing other numerical methods. Comparisons of 
solutions for linear static problems involving axisymmetric and unsymmetric 
loading, as well as nonlinear, axisymmetric problems, are availabe in the 
literature [99 , 102 ^ and will not be presented herein. 

The first set of problems to be investigated are static stability prob- 
lems which will be analyzed by employing the solution technique presented in 
Chapter k. This technique was first applied to problems involving cylinders. 
It was established that the technique produced accurate and rapid results 
using coarse structural idealizations [l57 ] • Difficulty was not encountered 
in predicting load reversals, or obtaining higher eigenvectors^ or eigenvectors 
which contained many waves within a segment. Thus, it is apparent that it is 
the number of integration points in a numerical integration method that is 
significant, rather than the number of segments. However, compared to finite 
element or finite difference methods, the segmentation utilized in the 
numerical integration method permits the use of much smaller matrices in the 
eigenvalue problem. 

Numerical Examples - Problem 1 : The first test problem involved the short 

cylinder shown in Figure 2k, subjected to an axial compressive end-load. The 

problem was chosen using the sizing parameters of Reference 142 to insure a 

critical mode shape in which the cross section remains circular (n=0), and 

there is only one half-wave along the length of the cylinder (m=l). The 
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Figure 2k Short End-Loaded Cylinder 
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boundary conditions used were classical simple supports, i.e., radial 
deformation was unrestrained until the point of incipient buckling. Tables 1 
through 3 show different aspects of the results. 

The first table demonstrates the accuracy of the results and the speed of 
convergence. As can be seen, only two iterations are necessary for the 
prediction of the first root. The second table demonstrates the speed of 
convergence to higher eigenvalues due to the additional information obtained in 
the current method. For a first trial load value, approximations to higher 
eigenvalues are available as well as to the lowest. As can be seen from the 
table when the first eigenvalue is obtained, a good approximation is available 
for the second eigenvalue and it can be obtained with only one additional 
iteration. Although this capability may not be overly important in stability 
analysis, it is very useful in free vibration analysis. The third table 
demonstrates the various capabilities of the current method with a coarse grid. 
As can be seen, no difficulty was encountered in obtaining eigenvalues corres- 
ponding to eigenvectors with many waves within a segment. Only two iterations 
were used to obtain each value in Table 3, and therefore the values should not 
be considered as fully converged. The first entry in the table was used to test 
the sensitivity of the method to negative eigenvalues. 

The last entry in Table 3 shows that an eigenvector with 15 half-waves in 
one segment was calculated correctly. The segment, of course, represents 8 
degrees of freedom in the stiffness matrix. To calculate such an eigenvector 
correctly would require up to 4 nodes per half-wave in a finite element 
idealization. Thus the equivalent degree of freedom requirement would be of the 
order of 4x15x4+4 = 244 d.o.f. Even in numerical integration the number of 
integration points must be kept to a reasonable limit for time considerations. 
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Table 1: Buckling of Short Axially Loaded Cylinder 



) 

Current Method (n=0) 



2 Segments 

1+ Segments 


Trial Load 

Result 

Trial Load ' 

Result 


1 x 10 ** 

32.1+5 x 10 1 * 

1 x 10** 

32.11+ x 10 ** 

32.07 x 10** 

32 x 10** 

32.09 x 10** 

32 x 10** 

32.09 x 10** 


Table 2: Short Cylinder Buckling Load Convergence 


Timoshenko 

(Ref.ll+2) 

Current Method (n=0) 1+ Segments 

1st 

Root 





Cumulative 
No. of 
Iterations 


dna x\oo*c 

jra noox, 


Trial 

Result 

Trial 

Result 

Trial 

Result 


1x10** 

32.11+xlO** 

1x10** 

88 . 92 x 10 ** 

1x10** 

200 . 18 x 10 ** 

1 

32.07x10** 

32x10** 

32.09x10** 

32x10** 

88.68x10** 

32x10** 

199.22x10** 

2 

88.12x10** 



88x10** . 

88.2x10** 

88x10** 

194.12x10** 

3 

193.1+3x10** 





191+xlO 1 * 

193.6x10** 

1+ 



































Table 3: Short Cylinder, High Buckling Loads 


Timoshenko (Ref. 142) 

Current Method (n=0) 

2 Segments 


Prediction 

% Difference 

X, = 32.07 x 10 k 

32.09 x hA * 

.06 

\ b = 769.26 x kA 

771.0 x 10^ 

.23 

X 7 = 1046.86 X 10^ 

1050.1 x 10 1 * 

.31 

A 8 = 1367.21 x 10 k 

1370.45 x 10 U 

.24 

X 9 = 1730.29 x 10 1 * 

1750.8 x^ 10 4 

1.2 

X/f- ^806 x hA 

5010 x 10 k 

4.2 

X /8 = 6920 x 10 H 

7321 x 10 k 

5.8 

\ zo = 8544 x 10* 

8680 x IQ 1 * 

1.6 

X 1S ~ 13,350 x 10^ 

14,410 X 10^ 

7.9 

X zg = 16,746 x 10 u 

17 ,212 x 10 k 

2.8 

\ 30 = 19 ,224 x 10^ 

19 ,080 x 10^ 

.75 


* Starting trial value was set at 1 x 10^ tension. 
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From the few test problems for this purpose it was noted that about 10 
integration points is a conservatively sufficient number to accurately 
describe one half-wave in a prospective eigenvector. 

Numerical Examples - Problem 2 : The second test problem involved a large, 

ring-stringer eccentrically reinforced cylinder (see Figure 25). The loading 
was a fixed internal stabilizing pressure of 31 psi. , in combination with a 
variable end load. Classical simple support boundary conditions were again 
utilized to enable comparison with References l4l and I50. The idealization 
used consisted of 20 segments for the whole structure, and 4 segments for the 
panel. Comparisons with analytical results for the overall and panel buckling 
modes are presented in Tables 4, 5 and 6. The overall critical mode was found 
to be n=0, m=13. Table 4 shows the analytical results for the n=0 calculations, 
and it can be seen that for this problem also, the convergence characteristics 
are excellent. By the fourth pass, the change from anticipated to corrected 
.value of the critical load is only .00017$. A comparison of the STARS-2B 
answers from the converged (fourth) pass, for estimates of some of the higher 
loads, shows an average difference of .57$ (with a maximum of 2.0$) for the 
first 11 roots, and an average difference of 2.72$ (with a maximum of 10$) for 
the first 18 roots, when compared to HASA THD 2960 (Ref. 150). Thus, when the 
first root is converged, excellent estimates are available for a large number 
of the higher roots. Similar results were obtained for the panel modes as 
shown in Table 5 • In this case convergence is obtained In two passes, although 
this is not verified until the third pass. The answers from References l4l and 
150 and the current work, should be close, but do not have to agree exactly, 
due to certain theoretical differences in the formulations. Reference 150 
uses Ponr.ell shell theory, while the current effort utilizes a Love-Reissner- 
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Table 4: Overall Buckling of Reinforced Cylinder 


1 


Root 

No. 

m 

NASA 
TND 2960 

NASA 

Current Method ( 20 segments ) 

CR 1280 

gnmn 

Trial 2 ( 6017.0) 

Trial 3 ( 5841.0) 

^^1^(5848.0) 

l=crit . 

13 

5848.15 

5842.9 

6017.77 

5840 .9 

5848.3 

5848.01 

2 

in 

5877.97 

5872.2 

6073.60 

5871.1 

5879.6 

5879.22 

3 

12 

5974 .05 

5968.4 

6126.76 

5972.4 

5978.8 

5978.56 

4 

15 

6032.72 

6025.6 

6264.54 

6033.4 

6043.1 

6042.75 

5 

16 

6290.95 

6281.6 

6445.00 

6308.4 

6315.4 

6315.67 

6 

n 

6301.05 

6293.9 

6569.61 

6310.3 

6319.6 

6319.13 

7 

17 

6637.3 

6624.8 

6971.53 

6683 *3 

6695.9 

6695.42 

8 

10 

6898.31 

6888.4 

7040 .98 

6923.9 

6928.9 

6928.68 

9 

18 

7060.56 

"7044.1 

7448.06 

7148.8 

7162.6 

7162.07 

10 

19 

755.2.39 

7530.9 

7940.50 

7691.5 

7704.7 

7704.21 

n 

9 

7875.44 

7861.1 

8023.79 

7923.9 

7928 .1 

7927.93 

12 

20 

8106.50 

8078.8 1 

9574.56 

8639.5 

8683.2 

8681.50 

13 

21 

8718 .03 


9821.53 

j 

9180.8 

9214.4 

9213.05 

l4 

22 

9383.22 


10196.41 

9492 .9 

9496.4 

9496.24 

15 

8 

9414.65 

9393.2 

11002.30 

10005.7 

10039-9 

10038.53 

16 

23 

10099.1 


12008.57 

10949.3 

10986.6 

10985.08 

17 

2k 

10863.3 


12010.31 

11950.2 

11952.7 

11952.57 

18 

25 

11674.00 


13167.73 

12009 .9 

12051.4 

12049.73 


n=0 = number of circumferential waves 
m = number of longitudinal half' -waves 
Results are single precision IBM 360/75 
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Table 5 : Panel Buckling of Reinforced Cylinder 


Root No. 

in 

NASA 
TND 2960 

NASA 
CR 1280 

Current Method (4 segments ) 

Trial-^lxlO^) 

Trial2( 51+80.0) 

Trial 3 ( 51+71.0) 

l=critical 

i 

5^93.93 

5I+I+9.6 

51+80 .7 

5470.3 

5470.3 

2 


6201.83 

6182.9 

6325.5 

6280.9 

6281.0 

3 

1 

10965 .00 

10881.3 

11314.1+ 

11168.0 

11168.1+ 

It 

B 

18626.90 

18372.1 

22557.5 

21631.9 

21633.7 


n=10 - number of circumferential waves 
m = number of longitudinal half-waves 
Results are single precision IBM 360/75 


Table 6: Buckling of Reinforced Cylinder, Reduction Scheme 


Current Method (Guy an reduction 20 segments >» b regions) . 

Root No. 

m 

Triali,.(581+8.0) no reduction 

Trialj t (581+8.0) reduction 

% difference 

1 

13 

5848.01 

581+8.01 

0.0 

2 

14 

5879.22 

5957.65 

1.3 

3 

12 

5978.56 

6280.32 

5.0 

1+ 

15 

60 1+2. 75 

65I+6.69 

8.3 

5 

16 

6315.67 

7065.60 

11.9 

6 

11 

6319.13 

7302.11 

15.5 

7 

IT 

669 5.42 

792 5. 61+ 

18.1+ 

8 

10 

6928 .68 

8130.09 

17.3 
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Kempner accuracy shell theory. Reference l4l on the other hand, while utilizing 
basic Love-Reissner theory, does not simplify some Fliigge accuracy terms such as 
(l ± £/R) when applying compatibility between rings or stringers, and the base 
shell. In addition, out-of -plane bending and twisting terms of the stiffeners 
are included therein. 

It should be noted that this problem contains closely spaced eigenvalues 
(see Table k). In the search of harmonics n = 1 and n = 2 the eigenvalues were 

t'M 

even more closely spaced. The current method encountered no difficulties in 
any of these cases. 

The overall shell buckling problem was also run using the matrix 
reduction scheme currently in the STARS -2B, -2V programs. The results are 
shown in Table 6. As can be seen, reducing the 20 segments to 4 regions does 
not affect the lowest root predictions, but does decrease the accuracy of the 
estimates for the higher roots. Thus, a consistent Guyan scheme can be used to 
analyze problems where idealizations demand a large number of segments. It is 
recognized that this scheme is basically for the reduction of static stiffness 
matrices, and other reduction schemes (for modal-eigenvalue problems) should be 
studied. However, the results of the current test problem serve to show the 
applicability of even the simpler reduction scheme within the accuracy of the 
STARS framework. 

It is interesting to qualitatively compare the above convergence' 
characteristics with those of B0S0R3 for a similar, stiffened cylinder problem 
(Ref. 138). In the STARS-2B stiffened cylinder problem, the analysis was 
started with an overall buckling load estimate of 17 . 1 $ of the converged 
critical load, and in four iterations the successive guesses were within .00017$ 
of each other. In the STARS-2B panel problem, the corresponding numbers were 
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18.28$, and in three iterations, results were within .0128$ of each other. In 
a similar stiffened cylinder problem solved with B0S0R3, using 91 finite differ- 
ence stations, and starting with a buckling load estimate of 99-5$ of the con- 
verged critical load, there is no convergence in single precision on the Univac 
1108 computer. Use of double precision produced convergence, and after five 
iterations the successive guesses are within .00383$ of' each other. It must 
also be noted in the comparison that the STARS -2B technique also provided 
excellent estimates to a large number of higher roots, while B0S0R3 found only 
the single lowest critical value. The difference in accuracy and speed of 
convergence, as well as the results provided (single or many roots) by each of 
the methods, is due to two factors. The major difference in the quality of the 
results is the fact that the matrices generated by the current numerical 
integration procedure are more accurate than those obtained by either finite 
differences or the finite element method. The number of roots immediately 
available is simply the result of using different numerical eigenvalue solution 
schemes. 

Numerical Examples - Problem 3 - As shown in Figure 26, the third stability 
problem studied in the present investigation is the PS-9 prolate spheroid, 
tested experimentally at the David Taylor Model Basin [153 ] * For this shell, 
a variety of theoretical results are available, and are tabulated in Fig. 26. 

The solution of Mushtari and Galimov [wd is based on the assumption that many 
lobes develop in both the circumferential and meridional directions. As 
apparent from the experimental results (n = 3), this is evidently not the case 
for this shell. The error in the theoretical predictions of Reference 153 is 
probably due to the assumption that the buckling deformation is confined to a 
narrow equatorial band of the shell. 
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Present Investigation 

210.6 

157.1*, 138.39+, 139.89'f , 135*1***- 

173.8 

DTMB Experimental Results (Ref. 153) 

- 

- 

DTMB Theory (Ref. 153) 

> 197 

197 

197 

Cohen [117] 

208.8 

(139.3) 138.7 

174.0 

Kalnins [ll8] 

Mushtari & Galimov [154] 


139.23 

95*5 psi (no harmonic prediction) 



* only membrane prestress terms included 

t membrane prestress and live pressure field terms included 
| predeformation neglected 

** all consistent nonlinear terms retained 

*** all consistent nonlinear terms retained (prestress matrices calculated with 
double precision arithmetic) 

\ 

Figure 2 6 Hydrostatically Leaded Prolate Spheroid 

162 









In the present investigation, the effects upon the critical buckling load 

of various nonlinear terms in the Equations (4-6, 7) were studied. The first 

number in Fig. 26 (*) is the buckling load based on the assumption that 

only the membrane prestress terms (Ngp, N p) are significant. The 

second value (+) is the buckling load based on the inclusion of the pressure 

rotation terms (f 0B , f^g, f^ B in Eq, (4-3)) but with the assumption that 

ei n ) = = 0. The third buckling load ($) was calculated by retaining all 

w o ^o 

nonlinear terms involving pressure or prestress, and neglecting only initial 
deformation. The final values of the buckling load, (**) and (***), are based 
on retaining all terms in Equations (b- 6 , j) , where the (***) result also 
includes the effect of double precision arithmetic in the calculation of the 
stiffness matrices. It may be observed that the greatest effect is obtained 
from the inclusion of the pressure rotation terms, and that the other effects 
are negligible by comparison. This is not surprising in the present problem 
since all the load is in the form of pressure, and predeformation (rotation) 
is expected to be minimal. 

The buckling loads predicted in this investigation and those of Cohen [117] 
and Kalnins £ll8 J are based on numerical integration. They are in excellent 
agreement with the experimental results. It is therefore apparent, inasmuchas 
the classical buckling load is obtained by experiment, that this prolate 
spheroid is not imperfection sensitive. It was expected that the predicted 
buckling loads will not be identical due to different assumptions in the 
theoretical formulations. Cohen uses Novozhilov [71] shell theory. Notice, 
that the number in parentheses given in Fig. 26, is obtained by the use of a 
nonlinear prebuckling state. In the present investigation and that of Kalnins 
[118] the Love-Reissner [161 shell theory is employed. However, Kalnins 
neglected some nonlinear terms, whereas in this investigation all consistent 
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terms (lOO, 137 J are retained. 

Numerical Examples - Problem 4 : The last test problem involves the study of the 

effects of axisymmetric prestress upon the vibration characteristis of the 
spherical cap shown in Figure 27. The present analysis is compared to other 
available results in Figure 28. As can be see^ the present STARS-2V results 
agree well with the Ebner [169] calculations using the VALORS £28! program, 
except for the n = 0 harmonic where a substantial disagreement is found. 

Although Ebner claims qualitative agreement with Reference 170 wherein results 
are available for the free vibration of spherical caps with R/h = 100 and an 
18° half opening angle, the following comparisons show otherwise: 



Ebner 

STARS -2V 

Ref. 170 

Harmonic 

II 

1.49&I+ 

1.2135 

1.26^9 

0 

o D 

II 

1. 1730 

1.15775 

1.1907 

1 


1. 3^25 

1.3325 

1.4398 

2 

§ 

o D 

il 

1.6522 

1.6359 

1.8166 

3 

o 3 

li 

2.0629 

2.03^2 

- 

4 

As can be 

seen above. 

Ebner' s zeroth harmonic frequency is greater than that of 


Ref. 170, whereas the frequencies for all the other harmonics are smaller. 
Similarly Ebner 's C0 Q is greater than Ch,. Neither of the above two items are 
consistent with the present analysis or Ref. 170. A further analysis with the 
Cohen program [l66] has confirmed the STARS-2V results tabulated above. 

The discrepancy in the Ebner calculations may possibly be explained by 
erroneous boundary conditions at the apex. Setting all displacements equal to 
zero as well as the meridional rotation is satisfactory only for n ^ 2, at the 
apex. With this boundary condition utilized for n = 0, the STARS-2V program 
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2 

yields a value of = 2.29136. Thus for this problem the apex boundary 

condition has a substantial effect upon the frequency results. A similar error 
■was found in the Kalnins £ll8] buckling analysis of problem 3. However, the 
effect for the buckling load in that problem proved negligible (see previous 
discussion and Fig. 26). > 
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Fig. 27 Shallow Spherical Cap 
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P/P 0 


Fig. 28 Variation of Natural Frequencies of Spherical Cap 
with Nonlinear Prestress 
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APPENDIX A 


RESULTANT STRESS-STRAIN RELATIONS FOR STIFFENED SHELLS 

One method of obtaining the necessary relationships between stress 
resultants and strains for eccentrically reinforced shells is based upon an 
"equivalent energy " approach. The energy of the composite system in 
terms of stress resultants and strains is equated to the energy of an equiv- 
alent orthotropic shell. 

The strain energy of a circumferential ring is given by [ 107] 
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where E ^ and G ^ are the modulus of elasticity and the shear modulus of 
the material from which the ring stiffeners are made, and is the tor- 

sional constant of the ring stiffeners. We will now distribute the strain 
energy of each ring uniformly over one half the panel spacing on each side 
of the ring. If the rings are spaced a distance S ^ apart, the total energy 
per panel length in the <f> direction may be written in the form 
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Assuming that the stiffeners are bonded to the shell and substituting Equation 
(1-13) for the total strain , , we obtain 
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This may be rewritten as 
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Similarly, the distributed strain energy per panel length in the 6 direction 
of n meridional stiffeners spaced a distance Sg apart (n Sg = 2 Trr o ) is 
given by 
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The strain energy of the unstiffened shell along a panel length is 
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This expression may be rewritten as 
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Combining like terms of Equations (A- 3) , (A-4) and (A- 5b) and 
using the two-dimensional stress -strain relations, we can obtain 
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where (see Figure A-l) and subscripts 0 and cp indicate coordinate directions, and 

A^ , Ag s cross-sectional area of the rings and the meridional 
stiffeners, respectively; 

Cr , Cg = eccentricity of the rings and the meridional stiffeners, 
respectively, measured from the reference surface of 
the shell (inwards positive); 

= moment of inertia of the rings and the meridional stiffeners. 


hi* *8 


respectively, about the basic shell centroidal axis; 
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Figure A-l Ring Stringer Reinforcement 
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J , J a torsional constant of the rings and the meriodional stiffeners, 
Jt\ b 

respectively; 

S_ , S c = spacing of the rings and the meridional stiffeners, 

R b 

respectively. 

Equations (A-6) are the relations between the stress resultants and 
the components of strain and curvature for the ring and stringer reinforced 
shell. They may be employed in lieu of Equations (1-18) in cases where the 
spacing of the ring and stringer reinforcement is such that the smearing 
technique yields valid results. Equations (A-6) can be rewritten in the 
following abreivated form. 
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Notice that inasmuchas the assumption was made in the Love-Reissner-Kempner 


first order theory that M 


<p0 


, from the last of Equations (A -7) it is 


implied that k^g = k ^ . From physical intuition, this relation may be 
approximate if the effects of the ring and meridional stiffeners are not identical, 
that is, the terms » an{ ^ Gg Jg/ Sg , in the torsional constant 

are not equal. 

Although Equations (A-6) and (A-7) were derived for a ring and 
stringer -stiffened shell, they could be extended to other cases by a suitable re- 
definition of the coefficients K„, and D^. . For example, they may be ex- 
tended to stiffened sandwich shells with equal or unequal face sheets, or to ring 
reinforced shells with corrugated skin in the meridional direction. Equations 

for the K.., C.. and D.. coefficients for the aforementioned cases are de- 
ll ij ij 

rived in References 32, 108. 

A more general form of the Equations (A-7) for layered media may 
also be obtained: 
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These equations are derived in Reference 32 . They are referred to a 

surface about which N^g is independent of k^g, and M^g is not dependent 

upon y g „ These equations may be employed in cases involving shells with: 
V 7 0 o 

a) homogeneous, sandwich, or multilayered skin reinforced by waffles 
at an arbitrary angle to the 8 , (p coordinate system , 

b) semi-sandwich shells (skin + corrugations) . 

The appropriate expressions for the , C.. and D. ^ coefficients for each 

case may be found in Reference 32 . 

Por a general isogria reinforcement Equations (A-9c, f) must be revised 


as follows: 
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APPENDIX E 


EQUATIONS FOR A DISCRETE RING 

The most consistent method of analyzing shells reinforced with widely spaced 
rings is to consider the shell segments and the rings as discrete structural 
members subjected to the given external loads and to the required interface 
conditions. The theory utilized is an expansion of that due to Cheney [l68] to 
include vibrations. The ring centroids and shear centers are allowed to be off- 
set, and the ring-shell connection can be eccentric. The necessary ring matrices 
for the programs are defined as: 

[M] Ring mass matrix (ring coordinates) 

{L^} Ring centrifugal acceleration load matrix (global coordinates) 

[k^] Ring prestressed stiffness matrix (ring coordinates) 

{L^} Ring thermal load matrix (ring coordinates) 

[T^ ] Transformation matrix to attached shell joint and shell global 
coordinate system 

The ring matrices are first converted to the shell joint and coordinate 
system, 

(£,} = J [T^f [t R ] {Aj} + ['2ro oJ -J £l A f [Ij,] (B-l) 

and then stacked appropriately for segment and region joints. The thermal load 
matrix is computed on the basis of a linear (radial) thermal distribution for the 
ring. 

The necessary matrices are presented in detail below: 
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where the necessary notation is defined in Figure B-l 
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APPENDIX C 


SPECIAL APEX CONDITIONS 

The system of Equations (1-27) through (1-29) , as derived in Chapter 
1, contains a singularity at r Q = 0 . In order to eliminate this singularity, 
and establish a suitable set of equations applicable at the apex we will make 
the terms in parentheses in Equations (l- 14c, d, f) and (l-12f) evaluated at - ■ 

(p = 0 equal to zero , and tnus apply L' Hospital' S rule. This may be accom- 
plished if the following conditions are satisfied at (p - 0 


u, e + v = v, 0 


U = w , e = 0 


°V, e " w e = w e, e + w cp = 0 


(C-l) 


The strain displacement relations at a smooth apex may then be written as 
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In obtaining the above relations. Equation (1-8) was tised to eliminate r 


o,(p * 


The equilibrium equations at the apex become: 



Q (p + Q 0, 0 = “ l N (p0 “e " N 0 U (p ], 0 [ N (p w 0 " N <p 0 ] 
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In a fashion similar to that of Reference 20 , the variables of 
Equations (C-2, 3 ) may be expanded in a Fourier series in the 0 direction, 
which involve only the part of the series (2-2) with the primed amplitudes. 
Thus from Equations (C-l) and (C-3) boundary conditions for the different 
harmonics are obtained, while from Equations (1-27 ) through (1-29), using 
L* Hospital* s rule and the established boundary conditions, the appropriate 
differential equations are obtained [ 6,8]. 

For the axisymmetric case (n= 0) the following conditions are ob- 
tained from Equations (C-l) through (C-3) 



Applying these relations and L* Hospital' s rule to the Equation system (1-27) 
through (1-29) , we obtain at the apex: 

194 



0 


0 


N 


( 0 ) 

"I 

’<P 9, £ 

r l 

( 0 ) 


JP’JL 

'i 


= 0 


U, 


( 0 ) 

(p 


( 0 ) 


w, 


<L 


*(°) 

r 

me = 

r i 


(0) (0) 

N„ F, 

je Jk 


i 


(l+2e^ 

V ’ 


M 


(0) 

S2LS2L 


= o 


(0) 


V, 


SL 


w 

r 


(0) 


1 


+ ^ K 22 ~ v 6cp K ll* 


N^°^(l-y J + N^, 

cp K V e((r t (ff 




)■ 

(0) (0) 

N e= \ 

(0) (0) 

M q = m 

0 (p 


For the first antisymmetric harmonic (n= 1) the following conditions 
are obtained from Equations (C-l) through (C-3) [ 20] . 
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Applying these relations, 
through (1-29) we obtain, 
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and L' Hopital* s rule to the Equations (1-27) 
at the apex: 
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For the second harmonic (ns 2), the following conditions are obtained from 
Equations (C-l) through (C-3) [ 20] , ; 
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Applying these relations , and L' Hopital' s rule, to the Equations (1-27) 
through (1-29), we obtain, at the apex: 
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Finally for the other harmonics (n > 3) the following conditions are obtained 
from Equations (C-l) through (C-3) [ 20 ] . 
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Applying these relations, and L* Hopital* s rule to the Equations (1-27) 


through (1-29) we obtain, at the apex; 
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' It should be noted that the nonlinear terms in these equations may be 
obtained from Equations (2. 9a) and (2. 9b) , It must be emphasized that the 
above equations are valid only for a smooth apex (sin cp = 0, cos (p = 1) . 

This method may be employed to obtain the appropriate equations for more 
complex apices. For a non-smooth apex {(p = <p ~ ■=£ 0) : 

sin (p - sin cp^ 0 
cos (p = cos (p Q 1 

Equations (C-5) through (C-ll) are utilized in the following manner 
in the numerical procedure discussed in Chapter 3. The apex boundary con- 
ditions from Equations (C-4, 6, 8, 10) are set in displacement form into the 
boundary condition matrix for the structure. The apex Equations (C-5, 7, 9, 11) 
are used only in the first step of the Runge-Kutta numerical integration pro- 
cedure, when it is applied to the edge r Q = 0 of the segment adjacent to the 
apex. In subsequent steps the sets of Equations (1-27) through (1-29) are 
employed. This procedure, however, is rather cumbersome, and it will be 
used only when the applied loading varies rapidly near the apex. In other 
cases, the apex boundary conditions (C-4, 6, 8, 10) may be satisfied at a circle 
of a very small but finite value of r Q . The results obtained on the basis of 
this approximation will be satisfactory at points away from the apex [ 6 , 8] . 
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APPENDIX D 


CONVERSION OF U.S. CUSTOMARY UNITS TO SI UNITS 

The International System of Units (SI) was adopted by the 
Eleventh General Conference on Weights and Measures in i960. Conversion 
factors for the units used in this report are given in the following 
table: 


Physical quantity 

U.S. Customary 
Unit 


SI Unit 
(**) 

Length 

in.. 

0.0254 

meters (m) 

Stress modulus 

ksi 

6.895 X 10 6 

newtons/meter^ 

(N/nr). 

Stress resultant 

ibf/in. 

175-1 

newtons/meter 

(N/m) 

Temperature change 

°F 

5/9 

Kelvin (K) 


* Multiply value given in U.S, Customary Unit by conversion factor to 
obtain equivalent value in SI Units. 

** Prefixes to indicate multiple of units are as follows: 


Prefix 

Multiple 

giga (G) 

10 9 

mega (M) 

10 6 

kilo (k) 

10 3 

deci (d) 

10- 1 

centi (c) 

10- 2 

milli (m) 

10 -3 
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